AFRL-MN-EG-TR-2001  -7023 


Development  of  a  Conservative  Integration  Scheme  for  the 
Equations  of  Hydrodynamics  Using  Unstructured  Grids 


Michael  L.  Stokes 


Mississippi  State  University 
Engineering  Research  Center 
Box  9627 

Mississippi  State,  MS  39762-9627 


GRANT  NO.  FO8630-98-1  -0004 

April  2001 

FINAL  REPORT  FOR  PERIOD  May  1998  -  September  2000 


DISTRIBUTION  A.  Approved  for  public  release;  distribution  unlimited. 


AIR  FORCE  RESEARCH  LABORATORY,  MUNITIONS  DIRECTORATE 

Air  Force  Materiel  Command  United  States  Air  Force  Eglin  Air  Force  Base 


20010813  033 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway.  Suite  1204,  Arlington,  VA  22202-4302, 

and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-01 88),  Washington,  DC  20503.  - - - — - 

1  AGENCY  USE  ONLY  (Leave  blank)  I  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

April  2001  Final  May  1 998  -  September  2000 


4.  TITLE  AND  SUBTITLE 


Development  of  a  Conservative  Integration  Scheme  for  the  Equations  of  Hydrodynamics 
Using  Unstructured  Grids 


6.  AUTHOR(S)  Michael  L.  Stokes 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Mississippi  State  University 
Engineering  Research  Center 
P.O.  Box  6810 

Mississippi  State,  MS  39762-9627 _ _ _ _ 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES)  (Program  Mgr  Name  &  Ph  #) 

USAF  Research  Laboratory 

Computational  Mechanics  Branch  (AFRUMNAC) 

101  W.  Eglin  Blvd,  Ste#337 
EglinAFB,  Florida  32542-6810 


5.  FUNDING  NUMBERS 
Grant  #:  F08630-98-1-1004 
JON:  25020727  PE:62602F 

PR:2502 
TA:07 
WU:27 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


10.  SPONSORING/MONITORING  AGENCY 
REPORT  NUMBER 


AFRL-MN-EG-TR-2001  -7023 


11.  SUPPLEMENTARY  NOTES 


Program  Manager:  Dr  Kirk  Vanden,  850-882-3124  (x3351] 


_ Availability  of  this  report  is  specified  on  the  verso  of  front  cover. _ 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT  12b-  DISTRIBUTION  CODE 

DISTRIBUTION  A.  Approved  for  public  release;  distribution  unlimited. 


13  ABSTRACT' 

In  this  report  we  describe  the  work  in  developing  Eulerian  upwind  one-dimensional  and  two-dimensional  numerical  models  for  shock 
impact  problems  on  unstructured  grids,  as  well  as  exact  solutions  of  one-dimensional  uniaxial  strain  problems  used  to  validate  the 
numerical  solutions.  This  report  contains  some  fundamental  theory  used  to  derive  the  working  equations,  as  well  as  examples  used 
to  illustrate  the  various  solutions. 


14.  subject  term  unstructured  hydrocode  finite-volume  elastic  plastic  impact  15-  number  of  pages 

81 

I  16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  OF 
REPORT 
UNCLASSIFIED 

NSN  7540-01-280-5500  ~~~ 


18.  SECURITY  CLASSIFICATION 

19.  SECURITY  CLASSIFICATION 

OF  THIS  PAGE 

OF  ABSTRACT 

UNCLASSIFIED 

UNCLASSIFIED 

20.  LIMITATION  OF  ABSTRACT 

SAR 


Standard  Form  298  (Rev.  2-89) 
Prescribed  by  ANSI  Std.  239-18 


NOTICE 


When  Government  drawings,  specifications,  or  other  data  are  used 
for  any  purpose  other  than  in  connection  with  a  definitely 
Government -related  procurement,  the  United  States  Government 
incurs  no  responsibility  or  any  obligation  whatsoever.  The  fact 
that  the  Government  may  have  formulated  or  in  any  way  supplied 
the  said  drawings,  specifications,  or  other  data,  is  not  to  be 
regarded  by  implication,  or  otherwise  in  any  manner  construed, 
as  licensing  the  holder,  or  any  other  person  or  corporation;  or 
as  conveying  any  rights  or  permission  to  manufacture,  use,  or 
sell  any  patented  invention  that  may  in  any  way  be  related 
thereto. 


This  technical  report  is  releasable  to  the  National  Technical 
Information  Services  (NTIS) .  At  NTIS  it  will  be  available  to 
the  general  public,  including  foreign  nations. 

This  technical  report  has  been  reviewed  and  is  approved  for 
publication . 


FREDERICK' A  DAVIS,  Technical  Director 
Assessment  and  Demonstrations  Division 


Dr  Kirk  Vanden 
Program  Manager 


If  your  address  has  changed,  if  you  wish  to  be  removed  from  our 
mailing  list,  or  if  your  organization  no  longer  employs  the 
addressee,  please  notify  AFRL/MNAC ,  Eglin  AFB  FL  32542-6810,  to 
help  us  maintain  a  current  mailing  list. 

Do  not  return  copies  of  this  report  unless  contractual 
obligations  or  notice  on  a  specific  document  requires  that  it  be 
returned. 


Contents 


1  Introduction  5 

2  One-Dimensional  Considerations  7 

2.1  Theoretical  Aspects  .  7 

2.1.1  Equations  of  Motion .  7 

2.1.2  Uniaxial  Strain .  9 

2.1.3  Two-wave  response  using  the  Mie-Grueneisen  Equation  of  State .  29 

2.1.4  Single- wave  plastic  response  using  the  Mie-Grueneisen  Equation  of  State .  32 

2.1.5  The  Equation  of  State .  32 

2.2  Numerical  Aspects .  34 

2.2.1  Integration  Scheme  for  the  Conservation  Equations .  34 

2.2.2  Flux  Evaluation  .  37 

2.2.3  Eigenvalue  Calculation .  39 

2.2.4  Integration  of  the  Constitutive  Equations .  40 

2.2.5  A  Turning  Algorithm  .  41 

2.3  Numerical  Results .  42 

2.3.1  Hydrostatic  Test  Case .  43 

2.3.2  Elastic  Solution  using  the  Experimental  Hugoniot .  43 

2.3.3  Elastic  Solution  using  the  Mie-Grueneisen  Equation  of  State  .  43 

2.3.4  Two- Wave  Solution  Using  the  Mie-Grueneisen  Equation  of  State .  44 

2.3.5  Single- Wave  Plastic  Solution  using  the  Mie-Grueneisen  Equation  of  State .  45 

3  Two-Dimensional  Considerations  47 

3.1  Conservation  Equations .  47 

3.2  Constitutive  Equations .  48 

3.3  Integration  of  the  Conservation  Equations .  48 

3.4  Integration  of  the  Constitutive  Equations .  50 

3.5  Results  for  Two-Dimensional  Plane  Strain .  50 

4  Conclusion  52 

5  Recommendations  53 

6  Additional  Figures  54 


1 


List  of  Figures 


2.1  Normal  Shock  Terminology .  13 

2.2  Single  shock  wave  impact  configuration .  14 

2.3  Single  shock  wave  impact  configuration  for  the  elastic  response .  18 

2.4  Figure  illustrating  the  two-wave  Plastic  Solution .  25 

2.5  Figure  illustrating  the  “Turning  Method”  .  46 

3.1  The  control  volume  using  a  dual  medium  grid .  49 

6.1  Pressure  for  hydrostatic  test  case .  55 

6.2  Density  for  hydrostatic  test  case .  56 

6.3  Energy  for  hydrostatic  test  case .  57 

6.4  Pressure  for  elastic  test  case  using  the  Experimental  Hugoniot .  58 

6.5  Density  for  elastic  test  case  using  the  Experimental  Hugoniot .  59 

6.6  Energy  for  elastic  test  case  using  the  Experimental  Hugoniot .  60 

6.7  Deviatoric  Stress  for  elastic  test  case  using  the  Experimental  Hugoniot .  61 

6.8  Pressure  for  elastic  test  case  using  the  Mie-Grueneisen  equation  of  state .  62 

6.9  Density  for  elastic  test  case  using  the  Mie-Grueneisen  equation  of  state .  63 

6.10  Energy  for  elastic  test  case  using  the  Mie-Grueneisen  equation  of  state .  64 

6.11  Deviatoric  stress  for  elastic  test  case  using  the  Mie-Grueneisen  equation  of  state .  65 

6.12  Pressure  for  two-wave  test  case  using  the  Mie-Grueneisen  equation  of  state .  66 

6.13  Density  for  two- wave  test  case  using  the  Mie-Grueneisen  equation  of  state .  67 

6.14  Energy  for  two- wave  test  case  using  the  Mie-Grueneisen  equation  of  state .  68 

6.15  Deviatoric  stress  for  two- wave  test  case  using  the  Mie-Grueneisen  equation  of  state .  69 

6.16  Pressure  for  single-wave  test  case  using  the  Mie-Grueneisen  equation  of  state .  70 

6.17  Density  for  single- wave  test  case  using  the  Mie-Grueneisen  equation  of  state .  71 

6.18  Energy  for  single-wave  test  case  using  the  Mie-Grueneisen  equation  of  state .  72 

6.19  Deviatoric  stress  for  single- wave  test  case  using  the  Mie-Grueneisen  equation  of  state .  73 

6.20  Pressure  for  impact  of  Copper  Bar  at  1000  m/s .  74 

6.21  Initial  Configuration  of  notched  bar .  75 

6.22  Deformed  Configuration  of  notched  bar  with  Pressure .  76 


2 


List  of  Tables 


2.1  Material  properties  used  in  the  example  problems .  17 

2.2  Results  of  the  numerical  procedure  for  the  shock  impact  problem  for  state  2  using  the  Mie- 

Grueneisen  Equation  of  State .  24 

2.3  Result  of  numerical  procedure  at  states  2  and  3  using  the  Mie-Grueneisen  equation  of  state 

for  the  two  wave  problem .  32 

2.4  Result  of  numerical  procedure  for  states  2  and  3  using  the  Mie-Grueneisen  equation  of  state 

for  the  single  plastic  wave  problem .  33 

2.5  Exact  vs.  Computed  Results  for  Hydrostatic  case(State  2) .  43 

2.6  Exact  vs.  Computed  Results  for  Elastic  case(State  2)  using  the  Experimental  Hugoniot.  ...  44 

2.7  Exact  vs.  Computed  Results  for  Elastic  case(State  2)  using  the  Mie-Grueneisen  Equation  of 

State .  44 

2.8  Exact  vs.  Computed  Results  for  the  Two-wave  Shock  using  the  Mie-Grueneison  Equation  of 

State .  44 

2.9  Exact  vs.  Computed  Results  for  the  Plastic  Single- Wave  response  using  the  Mie-Grueneisen 

Equation  of  State .  45 


3 


NOMENCLATURE 


c0  = 
e  = 
et  = 
A  <j>  = 
5<j>  = 


e  = 

f,fi  = 
T  = 
G  = 
To  = 

#o  = 
= 

Qij  — 

V  = 

Q,  = 
q  = 
p  = 
Po  = 


5(U)  - 
S*  - 
= 

*5Xy  = 
U  = 
t/5  = 

VJ  = 

V  = 
Vo  = 


constant  in  the  equation  of  state 

internal  energy  per  unit  mass 

total  internal  energy  per  unit  mass 

variation  of  (j>  in  space 

variation  of  (j)  in  time 

deformation  rate  tensor 

deviator  of  the  deformation  rate  tensor 

plastic  strain 

flux  vector 

flux  function 

bulk  Shear  Modulas 

plastic  material  constant 

y-intercept  of  the  plastic  yield  function 

slope  of  plastic  yield  function 

spin  tensor 

hydrostatic  pressure 

heat  flux 

velocity  vector(u,v,w) 
density 

pre-shock  density 

constant  in  the  equation  of  state 

deviator  stress  tensor 

stress  tensor 

equivalent  or  von  Mises  stress 
source  terms 

deviatoric  stress  in  x  direction 
deviatoric  stress  in  y  direction 
deviatoric  shear  stress  in  xy  plane 
vector  of  conserved  variables 
shock  speed 
tensor  velocity 
specific  density 
pre-shock  specific  density 
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Chapter  1 

Introduction 


The  US  military  has  known  that  certain  munitions  will  suffer  premature  ignition  if  the  internal  pressure  wave 
generated  from  a  mechanical  insult  is  of  sufficiently  large  magnitude.  Current  numerical  solutions  often  do 
not  provide  adequate  estimates  of  the  pressure  magnitude  due  to  the  application  of  artificial  viscosity  needed 
to  maintain  the  stability  of  the  numerical  scheme.  Other  inaccuracies  may  be  inadequate  resolution  of  grid 
points  in  regions  through  which  shocks  propagate,  and  inaccuracies  due  to  a  massive  skewing  of  the  grid 
subject  to  large  deformations.  A  similar  example  is  found  in  the  impact  of  high-speed  munitions  designed 
to  penetrate  concrete  bunkers  and  explode  at  depth  below  the  surface.  In  both  examples  it  is  necessary  to 
properly  capture  the  magnitude  of  the  pressure  wave  created  by  impact  induced  shocks. 

One  approach  to  accurately  capture  this  phenomena  is  to  cast  the  equations  of  motion  in  a  conservative 
Eulerian  form.  This  approach  has  several  advantages; 

•  Mass  continuity,  momentum,  and  energy  are  strictly  conserved.  Other  formulations  such  as  the  Ar¬ 
bitrary  Lagrangian  Eulerian  or  ALE  approach  maintain  mass,  but  the  conservation  of  energy  is  not 
strictly  conserved. 

•  The  Eulerian  form  can  be  cast  in  multiple  coordinate  frames  including  the  material  frame  where  the 
viewpoint  of  an  observer  is  coincident  with  the  particle  velocity  of  the  material,  or  from  a  reference 
frame  at  a  fixed  location,  or  from  an  arbitrary  reference  location  including  one  which  is  a  function  of 
time.  This  freedom  allows  a  great  deal  of  flexibility  in  dealing  with  the  quality  of  the  discretization 
during  extensive  deformation  of  the  solution. 

However,  despite  obvious  advantages  of  this  approach,  the  Eulerian  form  has  not  been  embraced  by  a 
majority  of  researchers.1.  The  classical  methods  of  integrating  these  equations  such  as  the  Lax-Wendroff 
class  of  schemes  are  highly  oscillatory  in  the  neighborhood  of  discontinuities  without  special  treatment. 
Fortunately,  much  work  has  been  done  in  the  Computational  Fluid  Dynamics  community  which  can  be  re¬ 
invented  for  the  equations  of  hydrodynamics.  One  example  is  the  notion  of  “upwind”  schemes  which  provide 
essentially  non-oscillatory  results  with  sharply  defined  shock  fronts.  The  current  research  extends  the  upwind 
approach  to  accommodate  grid  advection  during  the  integration  phase  of  the  solution  while  maintaining  at 
least  second  order  accuracy  in  both  time  and  space. 

This  document  describes  the  work  in  developing  one-dimensional  and  two-dimensional  numerical  models 
for  shock  impact  problems,  as  well  as  exact  solutions  of  one-dimensional  uniaxial  strain  problems  used  to 

^ee  reference  [1]  for  a  survey  on  these  techniques 
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validate  the  numerical  solutions.  It  contains  some  fundamental  theory  used  to  derive  the  working  equations, 
as  well  as  examples  used  to  illustrate  the  various  solutions.  The  numerical  examples  are  illustrated  in  Chapter 
6  using  screen  captures  of  the  working  programs  to  illustrate  the  behavior  of  the  numerical  codes. 
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Chapter  2 

One- Dimensional  Considerations 

This  chapter  develops  the  theory  for  exact  and  numerical  procedures  of  the  uniaxial  strain  problem  with 
examples  used  to  illustrate  the  concepts. 


2.1  Theoretical  Aspects 

The  contents  of  this  section  are  derived  from  a  number  of  existing  sources  including  [2], [3], [4],  and  [5],  many 
of  which  are  out  of  current  publication.  As  the  result,  some  of  the  developmental  material  is  presented  to 
provide  a  reasonably  complete  treatise  on  the  topic.  Much  of  the  theory  was  provided  and  explained  from 
unpublished  notes  authored  by  Dr.  Davy  Belk  of  Eglin  AFB.  The  remaining  material  was  developed  by  the 
author  as  necessary  to  complete  the  requirements  of  the  contract. 

2.1.1  Equations  of  Motion 

The  equations  of  motion  are  comprised  of  a  set  of  equations  which  provide  for  the  conservation  of  mass, 
momentum,  and  energy  along  with  a  set  of  non-linear  constitutive  equations  which  define  the  material 
derivatives  of  stress  and  strain.  The  equations  are  closed  with  an  appropriate  equation  of  state  relating 
pressure  to  internal  energy,  velocity,  stress,  and  strain.  The  fully  three-dimensional  form  of  the  equations 
will  be  presented,  then  reduced  to  one  dimension  using  the  assumption  of  uniaxial  strain  with  elastic/plastic 
deformation. 

Conservation  Equations 

The  conservation  law  equations  are  identical  for  both  solid  and  fluid  mediums.  The  only  variation  is  in  the 
form  of  the  stress  tensor  and  any  source  terms  that  may  be  present.  These  equations  in  divergence  form  are 


d(peT) 

dt 


§£  +  (p*>i)j  =  0 

(2.1) 

+  (pvivj  -  (Tji)}j  =  pfi 

(2.2) 

+  ( peTVj  -  VktTjk  +  Q j),j  =  pvkfk 

(2.3) 
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In  Cartesian  components,  the  stress  tensor  aij  can  be  represented  as  the  sum  of  the  deviator  stress 
and  the  hydrostatic  stress  components  as 


aij  —  sij  P&ij 

(2.4) 

where  the  pressure  p  is  defined  as  positive  in  compression  to  match  our  usual  idea  of  pressure  which  is  related 
to  the  stress  tensor  by 

&kk 

-p=  — 

(2.5) 

This  quantity,  the  negative  of  the  mean  normal  stress,  is  not  always  equivalent  to  thermodynamic  pressure, 
but  will  be  considered  as  such  in  this  context. 

Defining  Dij  as  the  classic  deformation  rate  tensor  and  D[-  as  the  deviator  of  the  deformation  rate  tensor 

Dij  —  2  +  v^i) 

(2.6) 

Dij  =  D{j  —  -Dkk6ij 

(2.7) 

The  Prandtl-Reuss  expressions  relating  deviator  stress  rate  (corrected  for  material  frame  indifference)  to 
plastic  strain  rate  are  given  by 

^  =  n„.„  -  .*(!„  +  2G  (d'v  - 

(2.8) 

1  -  ^ 

(2.9) 

where  0  is  defined  as 

_  k  .  \  k  ^  1  plastic  loading 

p  =  - 7777—  where,  <  _ 

l  _|_  H . Ce.)  I  k  ->  0  Otherwise 

3  G  \ 

(2.10) 

and  the  spin  tensor  Qij 

is  defined  as 

^ij  ~  9  (Vi,3  ~  V3,i) 

(2.11) 

The  set  of  conservation  laws  Eq.[2.1-  2.3]  and  the  constitutive  equations  Eqs.[2.8  -  2.9]  when  combined 
with  an  equation  of  state  (Section  2.1.5)  provide  a  complete  description  of  the  equations  of  motion  for  an 
elastic  or  plastic  medium. 
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2.1.2  Uniaxial  Strain 


To  reduce  the  constitutive  Eqs.[2.8][2.9]  for  the  case  of  uniaxial  strain,  it  will  be  assumed  that  no  shear 
stresses  can  exist  and  that  deformation  is  only  allowed  in  the  direction  of  motion.  For  this  case  the  stress 
tensor  is  not  one-dimensional,  but  rather  assumes  symmetry  in  the  transverse  direction.  For  this  case,  the 
stress  tensor  is  given  by 


mi  \  = 


ax  0  0 
0  ay  0 


0  0  Gy  J 

The  negative  of  the  mean  stress  from  Eq.2.5  is  given  by 

P=  +  2<T„) 

The  stress  deviator  components  are  calculated  by  solving  Eq.2.4  for  s* 


(2.12) 


(2.13) 


r  2 


i)  = 


Wx  Vy) 
0 
0 


J(<TX  Gy) 


0 


0  -|sx 


0 

2S 

0  0 

Similarly,  the  deviator  of  the  deformation  tensor  is  given  by 


0 

0 


-\(px-ay)  J 


(2.14a) 

(2.14b) 


m  = 


r  2  du 

3  dx 

o 

0 


The  equivalent  or  von  Mises  stress  a  is  defined  as 


0 

I  au 

'3  dx 

0 


0 

0 

i  du 
3  dx 


(2.15) 


-  -  (Z 
a  =  I  -Sijt 

Plastic  loading  occurs  when  the  equivalent  stress  reaches  the  yield  function  H'(e) 


iiSij  J  —  ^  |sx  |  —  Wx  &y  | 


(2.16) 


and 


>  0 


(2.17) 


(2.18) 
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If  the  two  conditions  specified  in  Eqs.[2.17]  [2.18]  are  satisfied,  then  k  in  Eq.[2.10]  is  1,  corresponding  to 
plastic  deformation,  otherwise  k  =  f3  =  0  and  the  strain  rate  given  by  Eq.[2.9]  is  0. 

With  the  previous  definitions  of  stress  and  strain,  the  conservation  and  the  constitutive  equations  for 
uniaxial  strain  reduce  to 


dp  djpu) 
dt  dx 

d(pu)  d(pu 2  +  P  -  Sx) 
dt  +  dx 

d(peT)  d[u(peT  +  P-  Sx)] 
dt  dx 

dSx 

dt 

dl 

dt 


0 


du 

dx 


(2.19a) 

(2.19b) 

(2.19c) 

(2.19d) 

(2.19e) 


Note  that  the  coupling  of  e  in  Eqs.[2.19a-2.19d]  is  due  to  the  fact  that  /3  (Eq.[2.10])  is  a  function  of 
H'(e)  which  influences  the  stress  rate  in  Eq.[2.19d].  However,  for  purely  plastic  deformation  where  Hf(e)  is 
a  constant,  /3  will  be  a  constant  and  Eq.[2.19e]  decouples  from  the  equation  set. 

Part  of  the  early  research  work  was  in  casting  the  advection  Eqs.[2.19d,2.19e]  in  a  form  where  the  entire 
equation  set  could  be  solved  as  set  of  partial  differential  equations.  However,  the  advection  equations  would 
have  to  be  manipulated  such  that  the  left-hand  side  was  in  a  divergence  form  with  all  remaining  terms 
treated  as  source  terms  on  the  right-hand  side.  The  advantage  of  this  scheme  would  be  a  higher  level  of 
coupling  between  the  five  equations  since  the  eigenvalues  of  the  coupled  system  would  have  the  contributions 
of  the  constitutive  equations.  There  is  however  no  unique  method  for  determining  which  terms  should  appear 
on  the  left-hand  side,  though  this  issue  has  been  a  point  of  research. 

As  following  [6],  consider  a  generalized  advection  equation  of  the  form 


4>t  +u(f>x=  0  (2.20) 

which  when  multiplied  by  p  and  combined  with  the  continuity  equation  [2.19a]  can  be  expressed  in  the 
divergence  form 


(fx!>)t  +  «)*  =  0  (2.21) 

Note  that  Eqs.[2.19d,2.19e]  can  be  put  in  this  form  for  the  special  cases  of  f3  =  1  and  j3  =  0,  respectively. 
The  well  known  jump  condition  for  Eq.[2.21]  is  given  by 


Pl<Pl{ul  -  Us)  =  Pr(Pr{uh  -  U8)  (2.22) 

where  U8  is  the  speed  of  the  discontinuity,  and  the  R  and  L  subscripts  refer  to  the  right  and  left  states, 
respectively.  If  the  jump  condition  for  continuity  given  by 


Pl{ul  ~  U8)  =  Pr{ur  -  Us) 


(2.23) 
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divides  Eq.[2.22],  then  <j>L  =  (j)R  indicating  that  (j>  would  be  continuous  across  the  shock.  Fedkiw[6]  concludes 
that  if  (j)  is  discontinuous  across  a  shock,  then  Eq.[2.21]  would  be  an  inappropriate  form  for  generalized 
advection  equations. 

While  the  development  of  the  argument  is  correct,  it  is  not  obvious  that  the  conclusion  is  definitive  for 
the  current  application.  It  should  be  obvious  that  his  argument  holds  only  in  the  case  where  the  additional 
equations  are  part  of  the  weak  solution  resulting  in  jumps  relations,  which  for  the  current  application  is  not 
applicable.  Additionally,  his  arguments  are  only  valid  for  equations  which  are  truly  in  a  divergence  form, 
which  again  does  not  match  the  current  case. 

Another  approach  to  determine  an  appropriate  divergence  form  of  the  constitutive  equations  was  to 
consider  the  eigen  properties  of  a  candidate  form  to  determine  if  the  equations  remained  hyperbolic.  Most 
of  the  variations  studied  had  singular  transformation  matrices,  resulting  in  an  incompatible  form.  One  form 
which  did  not  exhibit  this  property  was  the  form  given  by  Eq.[2.21].  However,  the  eigenvalues  for  this  form 
were  the  same  as  for  the  conservation  equations  alone,  thus  this  form  would  be  of  no  additional  advantage 
in  a  flux  based  algorithm. 

The  upshot  of  the  current  work  is  that  solving  Eqs.[2.19d,2.19e]  with  the  conservation  equations  as  a 
set  of  partial  differential  equations  is  not  recommended.  It  will  be  shown  that  the  exact  solution  of  the 
conserved  variables  is  dependent  only  on  the  initial  state  of  stress,  being  independent  of  the  derivatives  of 
Sx  or  6.  This  provides  a  natural  decoupling  of  the  constitutive  equations  which  should  be  mimicked  by  the 
numerical  procedure.  From  a  philosophical  standpoint,  attempts  to  pose  advection  equations  in  a  contrived 
conservative  form  could  create  unphysical  artifacts  in  the  solution  of  the  equations. 

Exact  Shock  Solutions 

The  “genuine  solution”  of  the  partial  differential  equation  in  the  general  form 


dp  dQjp)  _ 
dt  dx 


(2.24) 


is  strictly  only  valid  where  the  derivatives  are  defined.  Therefore,  a  discontinuous  solution  cannot  be  said  to 
be  genuine  when  p  or  Q(p)  are  discontinuous.  The  concept  of  weak  solutions  is  used  to  extend  the  solution 
set  to  include  discontinuities.  The  solution  set  is  represented  by 


-Us{p]  +  [Q]  =  o 


where  the  square  brackets  are  used  to  indicate  the  jump  across  the  discontinuity,  i.e. 


(2.25) 


[p]  =  P2~Pl 

[<5]  =  <?2  -  Qi 


and  Us  is  the  shock  velocity  given  by 


(2.26) 

(2.27) 


Us 


dx 

_ds_ 

dt 

ds 


dx 

dt 


(2.28) 
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The  general  solution  given  by  Eq.[2.25]  as  applied  to  the  set  of  conservation  equations  Eq.[2.19a-2.19c] 
without  heating  or  other  body  forces  is  given  by 


£ 

II 

.g. 

(2.29) 

Us[pv)  =  [pV2]  ~  [cr] 

(2.30) 

U,  p(e  +  iv2)  =  p{e  +  ^v2)v  -  [av] 

(2.31) 

These  equations  are  the  Rankine-Hugoniot  Equations  for  uniaxial  strain.  To  develop  a  general  solution 
framework  for  specific  problems,  a  special  case  is  derived  where  the  reference  frame  Us  =  0  resulting  in 


PlVi 

p\v\-o\ 
Pi(e  1  +  2Vi )Vl 


P2V2 

p2v\  -  cr2 

/>2(e2  +  V)V2  -(T2V2 


(2.32) 

(2.33) 

(2.34) 


where  V\  and  u2  are  velocities  relative  to  the  stationary  shock.  Substituting  Eq.[2.32]  into  Eq.[2.34]  results 
in 


e2-ei  =  -K2-u22)- 


£l  |  °2 

Pi  P2 


Using  Eq.[2.33]  to  eliminate  velocities  the  above  equation  can  be  written  as 


(2.35) 


e2-e1  =  i(u1+a2)(U2-U1)  (2.36) 

This  equation  is  often  referred  to  as  the  Hugoniot  equation.  Note  that  the  expression  is  independent  of 
coordinate  system  and  involves  only  thermodynamic  variables.  Another  useful  expression  is  found  by  solving 
for  U8  in  Eq.[2.29]  and  substituting  the  result  into  Eq.[2.30]  resulting  in 


(vi-v2)2  =  (*2-oi){V2-Vl)  (2.37) 

For  a  normal  shock  (Figure. 2.1)  propagating  with  unknown  velocity  U8  into  condition  2,  there  are  3 
Rankine-Hugoniot  equations  ,  but  there  are  five  unknowns  ,  i.e.  V2,P2,e2,  P2,  and  U8.  In  order  to  close  the 
system,  two  more  conditions  must  be  supplied.  The  conditions  are  typically  an  equation  of  state  relating 
pressure,  energy,  and  density,  together  with  any  one  of  the  state  variables  downstream  of  the  shock.  With 
this  information,  all  other  unknowns  at  state  2  can  be  determined. 

The  shock  problem  is  defined  as  the  impact  of  two  infinite  length  bars,  which  in  general,  may  be  different 
materials,  and  impact  one  another  at  a  contact  discontinuity.  Uniaxial  strain  will  only  be  considered,  thus 
the  bars  are  confined  in  such  a  way  as  to  prohibit  lateral  or  non-axial  deformation. 
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Figure  2.1:  Normal  Shock  Terminology. 


Experimental  Hugoniot 

Sometimes  a  Hugoniot  curve  is  erroneously  referred  to  as  an  “equation  of  state”  for  the  material.  An  equation 
of  state  is  an  entire  surface  in  (P,  V,e)  space  such  that  given  arbitrary  values  of  two  of  the  parameters,  the 
third  can  be  determined.  The  Hugoniot  is  a  single  curve  on  this  surface.  By  assuming  a  particular  form  of 
the  equation  of  state,  an  experimentally  determined  Hugoniot  is  often  used  to  determine  coefficient  values 
that  are  then  applied  off  the  Hugoniot.  The  misuse  of  the  “equation  of  state”  terminology  is  that  in  the 
absence  of  an  equation  of  state,  any  condition  involving  one  or  more  of  the  unknown  values  can  be  used  to 
make  the  system  determinate  as  long  as  the  condition  is  independent  of  the  jump  conditions.  This  condition 
is  then  used  as  a  replacement  for  the  equation  of  state  for  solving  shock  problems,  but  it  should  not  be 
mistaken  as  equivalent  to  an  equation  of  state. 

One  such  condition  is  a  linear  relation  between  shock  speed  Us  and  particle  velocity  v  given  by 


Us  =  Co  4*  sv  (2.38) 

where  the  reference  frame  is  such  that  the  initial  particle  velocity  is  zero  and  the  shock  speed  is  positive.  The 
coefficients  c$  and  s  are  determined  experimentally  where  cq  is  close  to  the  linear  longitudinal  wave  speed. 
For  the  case  as  shown  in  Figure  [2.1]  where  the  shock  speed  is  shown  positive  to  the  right,  then  Eq.[2.38]  is 
rewritten  as 


Us  ~  v2  =  Co  +  s(ui  -  v2)  (2.39) 

This  equation  is  called  the  right-running  wave.  If  the  wave  is  propagating  in  the  negative  direction  with 
respect  to  the  material,  then  the  formula  would  read 


Us  -v  1  =  Cq  +  s(v2  -  Vi) 


(2.40) 
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This  equation  is  called  the  left-running  wave.  One  of  these  empirical  relations  Eq. [2.39,2.40],  as  appropriate, 
can  be  used  to  analytically  represent  the  Hugoniot  curve  and  will  be  referred  to  as  the  Experimental  Hugoniot. 

Hydrostatic  Solution 

For  the  hydrostatic  case,  the  deviatoric  stress  is  set  to  zero.  The  impact  of  two  bars  will  initiate  a  single 
shock  wave  through  each  material  as  illustrated  in  Figure  [2.2].  The  left  material  will  contain  the  reflected 
shock  with  velocity  Urs  and  the  right  material  will  contain  the  transmitted  shock  with  velocity  Uts-  At  the 
impact  point  between  states  2  and  3,  a  contact  discontinuity  exists  with  velocity  Z7C,  therefore,  the  materials 
are  not  allowed  to  mix. 


Figure  2.2:  Single  shock  wave  impact  configuration. 

To  determine  which  properties  must  be  conserved,  consider  the  Rankine-Hugoniot  Equations  [2.29] 
through  [2.31]  as  applied  between  states  2  and  3  as  illustrated  in  Figure  [2.2].  Applying  Eq.[2.29]  results  in 


Us(p2-p3)  =  Uc(p2-p3)  (2.41) 

indicating  (1)  that  Us  =  UC)  or  simply  that  the  speed  of  the  discontinuity  must  equal  the  local  particle 
velocity  on  either  side  of  the  interface,  i.e.  Uc  =  v2  =  v 3,  and  (2)  that  a  jump  in  density  may  exist  at  the 
contact  discontinuity.  Applying  the  momentum  Eq.[2.30]  across  the  contact  discontinuity  results  in  a2  =  03, 
or  simply  that  the  normal  stress  must  be  continuous  across  the  contact  discontinuity.  For  the  hydrostatic 
case,  this  requires  p2  =  p3.  Applying  the  energy  Eq.[2.31]  across  the  interface  indicates  that  a  jump  in  energy 
may  exist  at  the  contact  discontinuity. 

For  the  configuration  shown  in  Fig. [2.2],  there  are  eleven  unknowns  including  Ur8i  v2,  p2,  e2,  and  p2 
from  the  left  material,  Uts,  u3,  p3,  e3,  and  p3  from  the  right  material,  and  Uc  from  the  contact  discontinuity. 
From  these,  three  are  duplicates  with  v2  =  i>3  =  Uc,  and  p2  =  p3,  therefore,  in  total  there  are  eight 
independent  unknowns.  For  the  left  material  there  are  three  Rankine-Hugoniot  jump  relations  and  one 
Experimental  Hugoniot  expression  relating  the  jump  conditions  between  states  1  and  2.  Similarly,  there  are 
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four  additional  conditions  from  the  right  material  relating  the  jump  conditions  between  states  3  and  4.  Since 
there  are  eight  equations  and  eight  unknowns,  the  problem  is  well-posed. 

The  solution  procedure  begins  by  finding  the  speed  of  the  contact  discontinuity.  This  provides  the 
velocity  at  states  2  and  3  by  satisfying  continuity  in  the  normal  stress  at  the  contact  discontinuity.  With 
this  information,  the  jump  conditions  will  be  analyzed  for  the  remainder  of  the  properties  at  states  2  and  3. 

The  density  at  state  2  is  expressed  using  Eq.[2.29]  as  a  function  of  the  respective  shock  speed  and  velocity 
jump  where 


92  =  Pi 


Urs  Vl 
Urs  ^2 


(2.42) 


and  similarly  for  state  3 


p3  =  p4 


Uts  ~  V4 

Uts  —  V3 


(2.43) 


Solving  the  momentum  equation  Eq.[2.30]  between  states  1  and  2  for  normal  stress  and  substituting  Eq.[2.42] 
results  in 


a2=cr1+  ( Urs  -  vi)(vi  -  v2)pi  (2.44) 

Similarly  between  states  3  and  4  we  have 

Oz  =  +  (Uts  ~  V4)(V4  -V3)p4  (2.45) 

Since  the  normal  stress  must  be  continuous  across  the  contact  discontinuity,  Eqs.[2.44,2.45]  are  equated  and 
simplified  by  substituting  expressions  for  the  shock  speeds  from  Eqs. [2.39,2.40]  to  give  a  quadratic  expression 
in  terms  of  Uc  where 


A2U2  +  A1UC  +  A0  =  0  (2.46) 

The  coefficients  of  the  polynomial  are  defined  as 

A2  =  sRp4  ~  sLpi 
A\  =  p\ (cq  +  2viSL)  4*  Pa(cr  —  2v4SR) 

Aq  =  p4 U4  (  CR  +  SRV4)  -  p\Vi  (Cq  +  SLV i)  +  cri  -  cr4 

where  the  superscripts  R  and  L  refer  to  the  right  and  left  materials.  For  A2  =  0, 

Uc  =  ~  (2.50) 

otherwise,  the  solution  is  given  by  the  negative  root  of  the  quadratic  formula 


(2.47) 

(2.48) 

(2.49) 
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(2.51) 


TT  —  Z—4  ~  V^-l  ~~  4AqA‘2 

2A2 

where  the  negative  root  is  required  such  at  v\  >  Uc  >  v4.  For  the  special  case  where  the  material  and 
thermodynamic  properties  are  identical  at  states  1  and  4,  Eq.[2.50]  reduces  to 

Uc=^(v  1+V4)  (2-52) 

The  shock  speeds  can  be  found  by  applying  Eqs.[2.40]  and  [2.39]  where 


Ure  =  Vl  ~  c£  +  SL(V2  ~  t>l)  (2-53) 

Uts  =  V4  +C„  +sr{v3  -  v4)  (2.54) 

The  densities  can  be  calculated  from  Eqs.[2.42]  and  [2.43],  and  the  pressures  using  Eqs.[2.44]  and  [2.45] 
realizing  that  p2  =  -02  and  p3  =  -er3.  Energies  can  be  calculated  from  Eq.[2.36]  where 

e2  =  e\  +  -(pi  +p2)(l/ 1  —  E2)  (2.55) 

e3  =  e4  +  i(p3+P4)(E4-E3)  (2.56) 

Note  that  this  solution  did  not  require  the  use  of  an  explicit  equation  of  state.  However,  it  can  be  shown 
that  the  Mie-Grueneisen  Equation  of  State  if  applied  to  the  hydrostatic  case  would  provide  the  same  pressure 
since  the  solution  stays  on  the  same  Hugoniot  line.  Also,  although  the  energy  equation  was  not  used,  it  can 
be  shown  that  the  solution  does  satisfy  the  conservation  of  energy  across  the  shock. 

EXAMPLE:  Hydrostatic  Response  using  the  Experimental  Hugoniot. 

Consider  two  bars  impacting  one  another  with  identical  material  properties  as  shown  in  Table  [2.1].  The  left 
bar  has  a  velocity  of  v\  —  100  m/s  and  the  right  bar  is  moving  at  v4  =  -100  m/s .  The  density  of  both  bars 
are  p\  =  p4  =  8930  The  pressure  and  energy  of  both  bars  are  zero.  The  problem  is  to  find  the  states  2 
and  3  as  shown  in  Fig.  [2.2]  assuming  a  hydrostatic  response. 

Using  Eq.[2.52]  to  calculate  the  velocity  of  the  contact  discontinuity 


Uc  =  0.5(100.0  +  -100.0)  =  0.0  m/s  (2.57) 

Using  Eqs.[2.53]  and  [2.54]  to  calculate  the  shock  speeds 

Urs  =  100.0  -  3940.0  +  1.49(0.0  -  100.0)  =  -3989.0  m/s 

Uts  =  -100.0  +  3940.0  +  1.49(0.0  -  -100.0)  =  3989.0  m/s 

Using  Eqs.[2.42]  and  [2.43]  used  to  calculate  the  densities 
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Table  2.1:  Material  properties  used  in  the  example  problems. 


Co 

= 

3940  m/s 

Po 

= 

8930  kg/ms 

r 

= 

2.0 

S 

= 

1.49 

H0 

= 

300  MPa 

H' 

= 

30  MPa 

E 

= 

121.9  GPa 

G 

= 

45  GPa 

P2 


P3 


8930 

8930 


(-3989.0  -  100.0) 
(-3989.0  -  0.0) 
(3989.0-  -100.0) 
(3989.00.0) 


9153.865634  kg/m3 
9153.865634  kg/m3 100 


The  normal  stresses  can  be  calculated  from  Eqs.[2.44]  and  [2.45]  as  follows 

<ra  =  0  +  (-3989.0  -  100.0)(100.0  -  0.0)  •  8930.0  =  -3.651477  GPa 
<r3  =  0  +  (3989.0  -  -100.0X-100.0  -  0.0)  •  8930.0  =  -3.651477  GPa 


Since  there  is  no  deviatoric  stress  for  the  hydrostatic  case,  the  pressures  at  states  2  and  3  are  simply  the 
negative  of  the  normal  stresses,  or 


P2  =  — <t2  =  3.651477  GPa 
P3  =  —<73  =  3.651477  MPa 

The  energies  are  calculated  using  Eqs. [2.55,2.56]  as  follows 

5000.0  J/Kg 
5000.0  J/Kg 


„  =  0.0  +  0.5.(0.0  +  3.651477e»)(^0-ssH5si)  = 
e5  =  0.0  +  0.5.(3.651477«»  +  0.o)(^0-s55^)  = 
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Elastic  Solution  using  the  Experimental  Hugoniot 

The  constants  c  and  s  given  for  the  Experimental  Hugoniot  will  not  be  applicable  for  impact  velocties  in  the 
elastic  range,  and  therefore  will  not  give  results  consistent  with  experiment1.  However,  it  is  instructive  to 
develop  these  relationships  for  the  reader  as  a  progression  of  complexity  both  from  the  mathematical  as  well 
as  science  perspective.  Hopefully  the  value  in  this  approach  will  be  appreciated  in  the  current  and  following 
sections. 

The  elastic  solution  is  similar  to  the  hydrostatic  solution  in  that  an  impact  will  generate  a  transmitted 
and  reflected  shock  as  shown  in  Figures  [2.3].  However,  unlike  the  hydrostatic  case,  deviatoric  stress  is 
produced  and  must  be  accounted  for  as  a  component  of  normal  stress. 


U  U  U 

rs  c  ts 


Reflected  Contact  Transmitted 

Shock  Discontinuity  Shock 


Figure  2.3:  Single  shock  wave  impact  configuration  for  the  elastic  response. 

There  are  fifteen  unknowns  which  are  Ur8,  p2,  v2,  P2 ,  ^2?  SX2,  and  <72  from  the  left  material  and  Uu , 
p3,  V3,  P3,  e3,  5X3,  and  for  the  right  material,  and  Uc  from  the  contact  discontinuity.  From  these  three 
are  duplicates  with  v2  =  v$  —  Uc  and  a2  —  cr3,  leaving  twelve  unknowns.  For  the  left  material  there  are 
three  Hugoniot-Rankine  Equations,  the  Experimental  Hugoniot,  a2  —  SX2  -p2y  and  Sx  =  Sx(p2 )  that  relates 
deviatoric  stress  to  density  resulting  in  six  independent  equations  for  the  left  material.  Similarly,  there 
are  six  equations  for  the  right  material  or  twelve  total  equations  matching  twelve  unknowns.  This  forms  a 
well-posed  system. 

The  calculation  of  the  contact  discontinuity  speed,  the  densities,  the  normal  stresses,  and  the  energies  at 
states  2  and  3  are  the  same  as  calculated  for  the  hydrostatic  case  (Section  2.1.2).  To  calculate  the  deviatoric 
stress,  the  assumption  is  made  that  the  material  response  is  elastic  and  that  the  Von  Mises  effective  stress 
<r  must  be  below  the  yield  surface,  or 

Hhe  use  of  the  Experimental  Hugoniot  will  severly  underpredict  the  shock  speeds  for  an  elastic  shock.  This  statement  also 
holds  for  the  hydrostatic  test  case  discussed  in  the  previous  section. 


(2.58) 


Since  the  value  of  deviatoric  stress  is  not  known  at  this  point,  the  procedure  is  to  simply  guess  that  the 
solution  is  elastic.  If  the  guess  is  wrong  then  the  elastic  procedure  will  be  abandoned  for  a  plastic  solution 
discussed  later  in  the  chapter. 

The  deviatoric  stress  is  determined  from  Eq.[2.19d]  cast  in  a  differential  form  where 

dSx  =  \Gde  (2.59) 

O 

where  for  an  elastic  response  0  =  0  and  =  e.  Substituting  the  natural  strain  formulation  de  =  j  in  the 
above  equation  and  integrating  from  some  initial  state  0  results  in 

/  ds-  =  lG)  7  <2'60) 

Sx  =  SXo  +  ^Gln  f-  (2.61) 

o  lo 

For  uniaxial  strain,  volumetric  and  longitudinal  strain  are  identical,  therefore, 


1  =  ^  =  ^ 

Jo  Vo  p 

Substituting  this  expression  into  Eq.[2.61]  leaves 

Sx  =  SXo  +  ^Gln  — 

o  p 

Writing  this  expression  between  states  1  and  2  shown  in  Fig. [2.2]  results  in 

SX2  =  SXl  +  t:Gl  In  —  (2.64) 

O  P2 

Similarly  for  state  3  the  deviatoric  stress  is  given  by 

5x3  =SX4  +  ^GR  In  (2.65) 

o  P3 

With  deviatoric  stress  calculated  and  validated  not  to  exceed  the  yield  surface  (Eq.[2.58]),  the  pressures  can 
be  calculated  at  states  2  and  3  by 


(2.62) 


(2.63) 
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P'2  =  Sx  2  -  (72 
P3  ^X3  &3 


(2.66) 

(2.67) 


The  energies  can  now  be  calculated  using  Eq.[2.55]  and  all  unknowns  are  defined. 

EXAMPLE:  Elastic  Response  using  the  Experimental  Hugoniot. 

The  same  configuration  for  the  material  properties  and  initial  conditions  will  be  used  as  in  Section  2.1.2, 
except  that  the  velocity  at  state  1  will  be  set  to  10  m/s  and  the  velocity  at  state  4  to  -10  m/s  so  that  the 
response  will  be  elastic.  Since  the  calculations  are  symmetric,  the  velocity  of  the  contact  discontinuity  Uc 
will  be  zero.  The  shock  speed  as  calculated  from  Eq.[2.53]  is 


Urs  =  10.0  -  3940.0  +  1.49(0.0  -  10.0)  =  -3944.9  m/s 
and  the  shock  speed  as  calculated  from  Eq.[2.54]  is 


Uts  =  -10.0  4-  3940.0  +  1.49(0.0  -  -10.0)  =  3944.9  m/s 
The  density  for  state  2  as  calculated  from  Eq.[2.42]  is 

*  =  8930U944489-o°o0)>  =  8952  637  ^ 

The  density  for  state  3  as  calculated  from  Eq.[2.43]  is 

”  -  8930<  (394499'-~ oT  =  8952-637  kS'm3 

The  normal  stress  at  state  2  as  given  by  Eq.[2.44]  is 


<r2  =  0.0  +  (-3944.9  -  10.0)(10.0  -  0.0)  •  8930.0  =  -353.1726  MPa 
The  normal  stress  at  state  3  as  given  by  Eq.[2.45]  is 

cr3  =  0.0  +  (3944.9  -  -10.0)(-10.0  -  0.0)  •  8930.0  -  -353.1726  MPa 
The  deviatoric  stress  at  state  2  is  given  by  Eqs.[2.64]  as 

s*2  =  0  +  l  ■  45.0  •  109  In  8g^°6°7  =  -151.9038  MPa 
and  the  deviatoric  stress  at  state  3  is  given  by  [2.65]  as 
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sxa  =  0  +  £  •  45.0  •  109  In  =  -151.9038  MPa 

3  3  8952.637 

The  energy  at  state  2  as  given  by  Eq.[2.36]  is 


62  =  »•»  +  0  5  '  (0  0  +  -353.1726*")  (55^5,  -  -^)  =  50.00  J,Kg 
The  pressure  at  state  2  is  given  by 

P2  =  SX2  —  02  —  —  151.9038e6  -  -353.1726e6  =  201.2688  MPa 
and  the  pressure  at  state  3  is  given  by 

p3  =  SX3  -  <r3  =  — 151.9038e6  -  -353.1726e6  =  201.2688  MPa 

With  the  deviatoric  stresses  known,  the  hypothesis  that  the  solution  is  elastic  can  be  tested  by  applying  the 
yield  criterion  Eq.[2.58] 


||5«|  =  ||  -  151.90381  =  227.8  MPa  <  a  =  300  MPa 

Since  the  effective  stress  is  below  the  yield  surface,  the  hypothesis  that  the  response  was  elastic  is  correct. 
All  other  properties  at  states  2  and  3  will  be  identical  to  the  hydrostatic  case. 

Elastic  Solution  using  the  Mie-Grueneisen  Equation  of  State 

This  second  method  of  calculating  the  elastic  response  utilizes  the  Mie-Grueneisen  Equation  of  State  which  is 
dependent  on  the  energy  equation.  The  problem  definition  is  the  same  as  in  the  previous  section  where  twelve 
unknowns  must  be  matched  with  twelve  independent  equations.  If  the  Mie-Grueneisen  Equation  of  State 
is  added  to  the  twelve  equations,  one  other  equation  must  be  removed  such  that  the  total  number  remain 
at  twelve.  The  only  possible  solution  is  to  replace  the  Experimental  Hugoniot  with  the  Mie-Grueneisen 
Equation  of  State. 

It  is  tempting  to  conclude  that  since  the  Mie-Grueneisen  Equation  of  State  derived  in  Section  [2.1.5]  uses 
the  the  Experimental  Hugoniot  in  its  derivation,  that  the  Experimental  Hugoniot  must  still  be  available 
for  the  calculation  of  the  shock  speeds  as  in  the  previous  example.  This  assumption  is  false.  It  is  true 
that  the  Mie-Grueneisen  Equation  of  State  as  implemented  calculates  reference  Hugoniot  values  using  the 
Experimental  Hugoniot.  However,  that  assumption  is  only  valid  within  the  context  of  the  derivation  of  the 
reference  Hugoniot  and  may  not  necessarily  be  available  outside  that  context. 

The  Mie-Grueneisen  Equation  of  State  has  a  very  complex  dependency  on  density  as  observed  in  Eq.[2.121]. 
This  complexity  makes  it  difficult  to  form  explicit  closed  form  solutions  for  all  but  very  simple  problems, 
thus  numerical  techniques  are  needed  to  algebraically  solve  for  the  unknown  states. 

To  calculate  the  jump  conditions  across  a  discontinuity  as  in  Fig.  [2.1],  consider  that  p2  and  all  parameters 
at  state  1  are  known.  From  the  Hugoniot  Equation  given  by  Eq.[2.36],  the  jump  in  energy  is  implicitly  solved 
giving 
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(2.68) 


ei  +  |  (ax  +  Sx(p2)  -  PH(V2)  [l  -  ^(V0  -  V2)])  (V2  -  V,) 

62  =  1  +  &(V2-Vt) 

where  Sx(p2)  is  given  by  Eq.[2.64]  and  P//(p2)  is  given  by  Eq.[2.120].  The  normal  stress  is  given  by 


a2  =  Sx(p2)  -  P(e2yV2) 

with  the  velocity  jump  solved  using  Eq.[2.37]  resulting  in 


(2.69) 


v2  =vx-  \J{p2  -  o\)(V2  -  V\) 


(2.70) 


where  the  negative  root  corresponds  to  a  left-running  wave.  For  a  right  running  wave 


v2=vx  +  y/(a2  -  ai){V2  -  Vi) 


(2.71) 


where  the  positive  root  corresponds  to  a  right-running  wave. 

The  difficulty  using  the  previous  analysis  is  that  specification  of  density  at  state  2  is  inconvenient  for 
solving  the  shock  impact  problem.  Ideally,  velocity  at  state  2  is  determined  by  satisfying  continuity  in  normal 
stress  at  the  contact  discontinuity,  then  using  the  velocity  to  solve  for  the  remaining  conditions  at  state  2 
and  3. 

To  accomplish  this,  a  numerical  procedure  is  defined  with  an  error  function  defined  as 


ffa)  =  v(p2)  -  Uc  (2.72) 

where  v(p2)  is  the  procedure  which  produces  velocity  in  Eq.[2.70]  given  density  p2  and  the  desired  contact 
discontinuity  velocity  Uc .  To  generate  a  new  density  which  reduces  the  error  in  Eq.[2.72]  from  some  initial 
estimate  of  density  p*,  the  following  function 


_  „  _  IM.  „  n  _  A  J  v(Pi)-Uc  \ 
Pi+1  ~  Pi  f'(Pi)  -  Pi  P{n(Pi  +  A p)  -  v(Pi)J 


(2.73) 


is  iterated  until  \pi+i~pi\  <  e.  The  method  is  functionally  a Newton-Raphson  technique  where  the  derivatives 
are  evaluated  using  a  first  order  finite  difference  approximation.  The  parameters  e  and  A p  are  reasonably 
small  values  taken  during  this  study  to  be  l.Oe"12  and  0.01,  respectively.  Convergence  usually  requires  two 
or  three  iterations  and  is  very  robust. 

To  calculate  the  contact  discontinuity  speed  Uc ,  Eq.[2.52]  is  used  if  the  material  and  thermodynamic 
properties  are  identical.  If  this  is  not  the  case,  then  the  result  of  the  previous  numerical  procedure  will  not 
produce  a  correct  normal  stress  at  the  contact  discontinuity,  i.  e.  a2  ^  <73  .  To  correct  this  problem,  a  second 
numerical  procedure  is  constructed  to  iterate  the  contact  discontinuity  speed  to  convergence.  Eq.[2.52]  is 
used  to  approximate  Uc  and  an  error  function  is  constructed  assuming  that  the  normal  stress  is  a  function 
of  velocity  where 
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f(v 2,^3)  =  ^3(^3)  -  cr2(v2) 


(2.74) 


With  this  function  a  similar  numerical  procedure  is  established  with  the  outside  loop  calculating  the  contact 
discontinuity  speed  and  the  inside  loops  calculating  state  properties  based  on  that  velocity.  However,  this 
choice  of  error  functions  is  not  tolerant  of  parameters  that  stray  far  from  the  symmetry  conditions.  An 
alternative  error  function  might  enforce  the  Hugoniot  Equation  Eq.[2.36]  at  the  contact  discontinuity  such 
that 


/(p2 5  P3)  =  e2(p2)  -  e3(p3)  -  -(cr2  +  cr3)(- - — ) 

2  P3  P2 


(2.75) 


though  this  function  has  not  been  investigated  thoroughly  to  report  on.  The  shock  speed  for  the  reflected 
shock  is  calculated  by  solving  Eq.[2.29]  for  Urs  resulting  in 


rr  _  pivi  -  P2V2 

urs  — 

Pi  -  P2 

Similarly  the  shock  speed  for  the  transmitted  shock  is  given  by 


Ut8  = 


P3V3  ~  P4V4 
p3  -  Pa 


(2.76) 


(2.77) 


EXAMPLE:  Elastic  response  using  the  Mie-Grueneisen  Equation  of  State. 

Consider  the  impact  problem  as  defined  in  the  previous  example.  The  contact  discontinuity  velocity  Uc  = 
0  m/s  will  be  the  same  since  the  configuration  is  symmetric.  The  results  of  the  numerical  procedure  as 
outlined  in  the  preceding  text  for  an  elastic  response  is  shown  in  Table  [2.2].  The  shock  speed  for  the 
reflected  shock  is  calculated  using  Eq.[2.76]  resulting  in 


Urs 


8930.0  •  10.0  -  8948.925  •  0 
8930.0  -  8948.925 


—4718.63  m/s 


Similarly  the  shock  speed  for  the  transmitted  shock  is  given  by  Eq.[2.77]  resulting  in 


Uts  = 


8948.925  •  0  -  8930.0  •  (-10.0) 
8948.925  -  8930.0 


=  4718.63  m/s 


Note  that  the  value  of  internal  energy  using  the  Mie-Gruneisen  equation  of  state  matches  that  for  the  elastic 
solution  given  by  the  Experimental  Hugoniot  even  though  the  pressure  and  density  vary  between  the  two 
solutions.  This  result  is  somewhat  unexpected  due  to  the  complexity  of  the  expression  for  energy  given  by 
Eq.[2.68]  as  compared  to  the  expression  using  the  Experimental  Hugoniot  given  by  Eq.[2.36].  This  suggests 
that  a  simplier  approach  might  be  available  for  the  solution  using  the  Mie-Gruneisen  equation  of  state,  though 
the  explanation  of  this  phenomena  remains  a  mystery.  This  point  certainly  needs  further  investigation. 

The  results  for  state  3  are  the  same  as  for  state  2  as  shown  in  Table  [2.2]. 
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Table  2.2:  Results  of  the  numerical  procedure  for  the  shock  impact  problem  for  state  2  using  the  Mie- 
Grueneisen  Equation  of  State. 


v-2  =  0  m/s 

p2  =  8948.93  Kg/m3 
p2  =  295.29  Mpa 
e2  =  50.0  J/Kg 
SX2  =  -126.98  Mpa 

aX2  =  -422.27  Mpa 


Two-wave  response  using  the  Experimental  Hugoniot 

When  the  effective  stress  exceeds  the  elastic  limit  on  the  yield  surface,  the  material  response  becomes  plastic 
and  is  accompanied  by  a  change  of  slope  on  the  stress-strain  curve  as  shown  in  Figure  [2.5].  When  this 
occurs,  two  shocks  may  form  depending  on  the  requirements  for  the  shock  speeds.  The  fastest  wave  is  called 
the  elastic  pre-cursor  and  compresses  the  material  to  the  elastic  limit.  The  second  shock,  if  it  exists,  is 
referred  to  as  the  “plastic  shock”  as  its  existence  is  due  to  the  plastic  response  of  the  material. 

For  this  case  the  impact  of  the  two  bars  results  in  six  states  as  illustrated  in  Figure.[2.4].  Though  it  is 
possible  for  one  material  to  have  an  elastic  response  and  the  other  a  plastic  response,  those  cases  will  not 
be  considered,  though  the  treatment  is  straightforward.  The  deviatoric  stress  generated  behind  the  elastic 
pre-cursor  at  states  1'  and  4'  reaches  a  maximum  defined  by  the  yield  surface  or 


2 -y  _  c  -  c 

'  Jl  —  <- ’x4l 


(2.78) 


where  the  negative  sign  assumes  compression.  Since  the  deviatoric  stress  is  known  at  states  1',  the  density 
can  be  calculated  from  Eqs.[2.64]  where 


Similarly  for  state  4' 


Pv  =  Pi  exp 


3 

4GL 


(SX1 


SXl,) 


(2.79) 


pv  =  pA  exp 


3 

4  GR 


{sX4-sXil) 


(2.80) 


The  velocity  at  state  1'  can  be  calculated  by  combining  the  expression  for  the  left-moving  Experimental 
Hugoniot  Eq.[2.40]  and  mass  continuity  Eq.[2.42]  to  provide 


V\ '  =  V\  + 


go jpv  -  Pi) 

sL(pv  -  Pi)  -  pv 


(2.81) 
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Reflected  Reflected  Contact  Transmitted  Transmitted 

Elastic  Plastic  Discontinuity  Plastic  Elastic 

Shock  Shock  Shock  Shock 


Figure  2.4:  Figure  illustrating  the  two-wave  Plastic  Solution. 


Similarly  for  state  4' 


V# 


zP±L 

SR(p4'  -  Pi)  -  Pi' 


The  shock  speed  Ure  is  calculated  from  Eqs.[2.40]  to  give 


Ure  =Vi  -Cg  +SL(Vy  -  Uj) 

Similarly,  the  shock  speed  Ute  is  calculated  from  Eq.[2.39]  to  give 


Ute  =  V4  +  Cq  +  sR(v 4,  -  v4) 
The  normal  stress  can  be  computed  from  Eq.[2.44]  at  state  1'  using 


(2.82) 


(2.83) 


(2.84) 


CTi-  =  <7i  +  pi  (Ure  ~  Vi)(vi  -  Vy) 


(2.85) 


and  similarly  at  state  4'  as 


cr4,  =  a4  +  Pa(Uu  -  v4)(v4  -  v4l) 


(2.86) 
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Pressure  at  state  V  is  calculated  by 


Pv  =  SXl,  -  av  (2.87) 

and  similarly  for  state  4' 

Pi'  =  Sx4,  -  a  A'  (2.88) 

Energy  can  be  computed  using  the  Hugoniot  Equation  [2.36]  resulting  in 

er  =e1  +  i(<71+<r1,)(Vi'-U)  (2.89) 

and  similarly  for  state  4'  where 

e4,  =  e4  +  +  (r4.)(V4.  -V4)  (2.90) 

Note  that  the  solution  for  state  1'  and  4'  does  not  depend  on  the  velocity  of  the  contact  discontinuity  as  in 
previous  examples.  The  states  at  1'  and  4'  are  completely  determined  by  the  elastic  yield  stress.  However, 
this  is  not  the  case  for  states  2  and  3. 

If  the  material  response  is  purely  plastic,  i.  e.  H '  =  0,  then  a  shock  is  necessary  because  of  a  thermo¬ 
dynamic  mismatch  between  states  1'  and  2,  and  similarly,  between  states  4'  and  3.  The  deviatoric  stress 
across  the  plastic  shock  will  be  continuous  since  the  bar  cannot  unload  and  the  yield  surface  constrains  the 
deviatoric  stress  from  increasing.  Given  the  speed  of  the  contact  discontinuity  £/c,  the  shock  speeds  Urp  and 
Utp  can  be  determined  from  the  Experimental  Hugoniot  relation,  and  the  rest  of  the  unknowns  at  states  2 
and  3  follow  the  same  procedure  as  for  a  hydrostatic  shock  (See  Section  [2.1.2])  where  states  1  and  4  in  the 
hydrostatic  solution  are  replaced  with  V  and  4',  respectively. 

For  a  material  with  Hf  >  0,  then  in  addition  to  the  thermodynamic  jump  described  in  the  previous 
paragraph,  there  would  be  an  additional  production  of  deviatoric  stress  from  plastic  compression  governed 
by  Eq.[2.19d].  The  yield  surface  itself  depends  on  the  production  of  plastic  strain  l  governed  by  Eq.[2.19e]. 
These  two  production  equations  can  be  combined  to  relate  stress  and  plastic  strain  similar  to  the  elastic 
case.  The  ratio  of  Eq.[2.19d]  and  Eq.[2.19e]  gives 


4S_x_. 

dt 

dt 

dt 


| <7(1  -  0)vx 


=  2 G(i  -  1)/  = 


dSx 

de 


(2.91) 


where  /  =  Sign(vx).  In  a  compressive  environment,  vx  <  0,  therefore,  /  =  — 1.  Rewriting  Eq.[2.91]  in  a 
differential  form  leaves 


dSx  =  2G(1  -  l)dc 


(2.92) 
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For  elastic  strain,  de  =  y ,  but  the  current  definition  of  elastic  strain  requires  that  e  >  0  which  is  violated  in 
compression.  To  compensate  for  this,  plastic  strain  will  be  defined  in  terms  of  l  where 


dl=-  j  (2.93) 

Substituting  Eq.[2.93]  into  Eq.[2.92]  leaves 

dSx  =  2G(^-l)j  (2.94) 

Integrating  Eq.[2.94]  from  the  elastic  yield  point  Sj  to  the  resting  point  on  the  yield  surface  provides 


/ 


sr 


dSx  =  2 G(i  - 


<ti 

I 


or, 


sx-s?  =  2G(i-l) 


1 1 


In  l 


-Wo 


Sx-S%  =  2G(I -lKlnf- Info)  =  2G(i-l) Ini 
As  in  the  elastic  case  for  uniaxial  strain,  volumetric  and  longitudinal  strain  are  identical,  therefore 


(2.95) 


i=  V  =po 
lo  Vo  P 

Substituting  the  equation  above  in  Eq.[2.95]  results  in 

S,  =  Sr0+2G(i-l)ln^ 

Rewriting  this  equation  in  the  context  of  state  2  from  Fig.  [2.4] 

SX2  =  SYx\L  +  2G{~-\)\nPj- 
and  for  state  3  the  equivalent  expression  is 


(2.96) 


(2.97) 


(2.98) 


However,  at  this  point  SX2  and  SX3  can  not  be  calculated  because  the  densities  at  those  states  are 
unknown.  However,  the  speed  of  the  contact  discontinuity  Uc  =  v 2  —  ^3  can  be  calculated  as  in  the  solution 
of  Eq.[2.46]  where  state  1  is  replace  with  V  and  state  4  with  4'.  The  shock  speeds  Urp  and  Utp  follow  suit 
from  Eq.[2.40]  and  Eq.[2.39],  respectively,  where 


Urp  =  VV  Cq  +  SL(V 2  -  Vv) 


(2.100) 


and 


Utp  =  v4'  -  c£  +  sR(v 3  -  v4')  (2.101) 

The  densities  at  states  2  and  3  can  be  calculated  from  Eq.[2.42]  and  Eq.[2.43]  as 


P  2 


=  pv 


Urp  vy 

U rp  ^2 


(2.102) 


and 


P3 


Utp  -  v4* 
_  PA'  Utp  -  V3 


(2.103) 


Now  that  the  densities  are  known  at  states  2  and  3,  the  deviatoric  stresses  in  Eq.[2.98]  and  Eq.[2.99]  can  be 
calculated. 

The  normal  stresses  at  states  2  and  3  are  calculated  from  Eq.[2.44]  and  Eq.[2.45]  as 


<J2  =  ay  +  pv(Urp  -  vv)(vy  -  V2) 


(2.104) 


and 


0-3  =  <?4'  +  Pa1  (Utp  -  V4')(v4>  -  V3) 


(2.105) 


Pressures  are  calculated  directly  as 


p2  =  SX2  -  <J2 


and 


P3  ^3 


and  energies  are  calculated  from  Eq.[2.55]  and  Eq.[2.56]  as 
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(2.106) 


e2  =  ev  +  -{pv  +P2){Vy  -  V2) 


and 


e3  =  e4<  +  i(p4'  +  PaWv  -  V3)  (2.107) 

An  example  for  this  case  will  not  be  given  because  while  it  has  academic  value,  a  two  wave  system 
given  the  current  material  properties  is  not  possible  without  an  unrealistically  high  yield  stress,  therefore 
the  example  will  be  left  to  the  interested  reader. 

2.1.3  Two- wave  response  using  the  Mie-Grueneisen  Equation  of  State. 

The  description  for  the  two- wave  response  is  very  similar  to  the  case  for  the  Experimental  Hugoniot,  except 
of  course  that  the  Experimental  Hugoniot  is  being  replaced  with  the  Mie-Grueneisen  Equation  of  state.  The 
deviatoric  stress  and  density  at  state  1'  and  4'  will  be  identical.  Once  the  density  is  known,  the  Hugoniot 
curve  Eq.[2.36]  is  expanded  in  terms  of  the  deviatoric  stress  and  pressure  at  state  1/.  Solving  implicitly  (See 
Eq.[2.115])  for  energy  at  state  1'  leaves  the  rather  lengthy  expression 


ev  = 


ei  +  1  (t,  +  SAP,')  -  ph(Vv )  [l  -  jftr(Vo  -  IV)])  'tv  -  Vi) 
i  +  sfcOV-H) 

The  pressure  at  state  1'  is  found  by  solving  the  Hugoniot  Equation  Eq.[2.36]  for  pressure  or 


(2.108) 


Pv 


+  (7\  *“ 


2  •  (ey  -  eQ 
Vv  ~  Vi 


The  velocity  is  found  from  Eq.[2.70]  as 


Vv  =Vi-  \/ (cri/  -  (Ti)(Vl-  -  Ki) 

where  ov  =  SXl,  -  pv  ■  The  velocity  for  state  is  found  from  Eq.[2.71]  as 


(2.109) 


(2.110) 


W  =V4  +  v/(<74<  -  <7 4)(^4'  -  VA)  (2.111) 

thereby  completing  the  unknowns  at  state  1'. 

The  states  at  2  and  3  must  be  determined  numerically  as  was  the  case  with  the  elastic  solution  using  the 
Mie-Grueneisen  Equation  of  State.  The  contact  discontinuity  speed  is  calculated  or  approximated  using 

Vc  =  \(v1+v4)  (2.112) 
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The  numerical  procedure  developed  in  Section  2.1.2  is  used  where  Sx(p2)  uses  the  plastic  function  as  defined 
in  Eq.[2.98].  The  states  at  2  and  3  are  calculated,  then  if  the  normal  stresses  are  not  equal  at  the  contact 
discontinuity,  the  densities  at  states  2  and  3  are  recalculated  using  a  Newton-Raphson  procedure  as  devel¬ 
oped  earlier  until  the  error  tolerance  is  below  the  established  criteria. 

EXAMPLE:  Two-wave  solution  using  the  Mie-Grueneisen  Equation  of  State. 

The  physical  orientation  will  be  the  same  as  in  previous  examples  and  the  material  values  will  be  the  same 
as  in  Table  2.1.  The  initial  properties  at  states  1  and  4  will  be  the  same  as  in  the  previous  examples  except 
that  v\  =  100  m/s.  The  problem  is  to  find  the  states  at  1',  2,  3,  and  4'. 

The  problem  is  the  collision  of  two  like  metal  bars  which  impact  on  another.  The  physical  properties  are 
given  in  Table  2.1  except  that  Hb  will  be  set  to  zero  such  that  the  response  is  purely  plastic.  The  initial 
conditions  in  the  bars  will  be  pi  —  Pa  =  0.0,  e\  =  e\  =  0.0,  SXl  =  SX4  =  0.0.  The  velocity  of  the  left  bar  will 
be  v\  —  100  m/s  and  the  velocity  of  the  right  bar  will  be  va  =  -100  m/s. 

The  assumption  is  made  that  the  deviatoric  stress  exceeded  the  yield  point  assuming  an  elastic  response, 
thus  the  solution  must  exhibit  a  plastic  response.  However,  to  have  a  two-wave  system  the  shock  speed  of 
the  plastic  wave  must  be  less  than  the  shock  speed  of  the  elastic  wave.  If  this  is  not  the  condition,  then  the 
shocks  will  coalesce  into  a  single  shock.  From  Eq.[2.78]  we  have  that  SXl,  =  SX4 ,  =  §  *  300  =  -200  Mpa  in 
compression.  From  Eq.[2.79],  we  can  calculate  the  density  at  state  1'  as 


pv  =  8930.0  •  exp 


3.0 


4.0  •  45. Oe9 


(0.0  - (-200.0 -106)) 


8959.816332  Kg/m 3 


The  calculation  of  energy  at  state  1'  will  be  broken  into  several  pieces  to  simplify  the  solution.  The 
reference  Hugoniot  ph(Pv)  is  given  by  Eq.[2.120]  or 

3940  02  •  ( — - - - - ) 

pH(pV)  =  -  18930  0 - 8959.81633a /  =  466.()77  MPa 

fl  49.  (— 1 _ I _ ) _ 1 — l2 

v  8930.0  8959.816332  /  8930.0  J 


The  energy  is  solved  using  Eq.[2.108] 


ev 


0.0  +  I  (o.o  +  -200  •  106  -  466.077  •  106 


1  - 


2.0 


& 8959.816332 


(-k _ l _ 

'  8930  8959.816332 


>])< 


1 


8959.816332  8930 


L— ) 

)30  f 


1  + 


2.0  / 
2  1  \ 
c 8959.816332 


1 


8959.816332  8930 


)30J 


ev  =  124.272  J/Kg 

The  pressure  is  found  at  state  V  using  Eq.[2.109]  where 


pv  =  —200.0  •  10  4-  0.0 


|6  .  nn  2.(124.272-0.0) 


i 


l 


:  466.744826  MPa 


8959.816332  8930.0 


The  normal  stress  at  state  1'  is 


ov  =  SXl,  -  pv  =  -200  •  106  -  466.744826  •  106  =  -666.744  MPa 
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The  velocity  is  found  at  state  V  using  Eq.[2.70] 


**  =  100.0  -  J  (  -  666.744e6  -  O.o)  ^ >  =  84-28  ^ 


The  shock  speed  Ure  is  calculated  using  Eq.[2.29]  where 


8930.0  •  100.0  -  8959.826  •  84.2347 
8930.0  -  8959.826 


—4635.95  m/s 


The  properties  at  state  4'  are  identical  to  state  1'  expect  the  velocity  which  is  calculated  from  Eq.[2.71] 
which  is 


vv  =  -100.0  + 


666.744e6  -  0.0 


1 


8959.826  8930 


—  )  =  —84.23  m/s 

i.O  J  ' 


The  shock  speed  Urp  is  calculated  from  Eq.[2.29]  where 


8959.826  •  84.2347  -  9147.919  •  0 
Urp  ~  8959.826  -  9147.919 

Summarizing  the  solution  at  state  1': 


-4012.535  m/s 


p  =  8959.816332  Kg/m3 
v  =  84.23  m/s 
p  =  466.744826  MPa 
Sx  =  -200  MPa 
a  =  —666.744  Mpa 
e  =  124.272  J/Kg 

and  summarizing  the  results  at  state  4' 

p  =  8959.816332  Kg/m3 
v  =  —84.23  m/s 

p  =  466.744826  MPa 
Sx  =  -200  MPa 
a  =  —666.744  Mpa 

e  =  124.272  J/Kg 

The  conditions  at  state  2  and  3  are  solved  using  the  described  numerical  procedure  with  the  results 
presented  in  Table  2.1.3. 
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Table  2.3:  Result  of  numerical  procedure  at  states  2  and  3  using  the  Mie-Grueneisen  equation  of  state  for 
the  two  wave  problem. 


p  = 

9147.919  Kg/m6 

V  = 

0  m/s 

p  = 

3.559  GPa 

sx  = 

-200.0  MPa 

a  = 

-3.7587  GPa 

e  = 

5219.884  J/Kg 

2.1,4  Single- wave  plastic  response  using  the  Mie-Grueneisen  Equation  of  State. 

If  the  contact  discontinuity  velocity  increases  in  a  two-wave  system,  at  some  point  the  velocity  of  the  plastic 
wave  will  equal  the  velocity  of  the  elastic  wave.  When  this  happens,  the  shocks  collapse  back  to  a  single 
wave  system  as  shown  in  Fig.[2.2].  The  production  of  deviatoric  stress  at  state  2  is  taken  from  Eq.[2.98] 
where 


SX2=S^\L+2G(^-l)ln^- 
P  P  2 


(2.113) 


and  similarly  for  state  3 


5*3  =  +  2G(-i  —  1)  In  —  (2.114) 

P  P3 

The  numerical  procedure  used  in  Section  2.1.2  is  used  to  generate  the  parameters  at  state  2  and  3. 

EXAMPLE:  Single  Wave  Plastic  Solution  using  the  Mie-Grueneisen  Equation  of  State. 

For  this  case  the  velocity  at  state  1  will  be  increased  to  v\  =  600  m/s  with  V4  —  -600.0  m/s.  The  numerical 
procedure  as  described  previously  will  be  employed  to  solve  for  the  unknown  states  at  2  and  3.  The  results 
of  this  procedure  are  shown  in  Table  2.4. 

2.1.5  The  Equation  of  State 

A  convenient  and  often  used  equation  of  state  for  shock  related  problems  is  the  Mie-Grueneisen  equation  of 
state,  which  gives  pressure  in  terms  of  internal  energy  and  specific  volume  or  density  as  follows: 

p(e,V)=pH(V)  +  ^(e-eH(V))  (2.115) 

Its  convenience  in  part  stems  from  the  use  of  the  shock  Hugoniot  to  give  reference  values.  The  pressure  at 
an  arbitrary  value  of  specific  volume  and  internal  energy  p(e,  V)  is  obtained  by  finding  the  pressure  Ph(V) 
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Table  2.4:  Result  of  numerical  procedure  for  states  2  and  3  using  the  Mie-Grueneisen  equation  of  state  for 
the  single  plastic  wave  problem. 


p  =  10189.95  Kgjn? 

v  —  0.0  m/s 

p  =  25.8  GPa 
e  =  .18045  MJ/Kg 

Sx  =  -200.0  MPa 

a  =  -26.0  GPa 


and  internal  energy  e//(V')  on  a  reference  Hugoniot  at  a  given  density.  The  variation  of  the  desired  energy 
from  the  Hugoniot  energy  is  then  linearly  related  to  the  pressure  difference  at  the  fixed  density,  with  the 
proportionality  constant  being  where  T  is  the  Grueneisen  coefficient. 

The  Hugoniot  reference  state  with  a  specified  density  pH  must  satisfy  the  Rankine-Hugoniot  relations 
with  a  “zero  state”  where  pressure,  energy,  and  velocity  are  zero,  i.  e.  po  =  e0  =  vo  =  0,  with  specified 
density  po.  Letting  1  state  in  Eq.[2.35]  be  the  “zero  state”  and  state  2  be  the  Hugoniot  reference  state,  then 


cn(Vff)  =  ^(V1 b-Vfr) 

Using  the  same  substitution  in  Eq.[2.37]  and  solved  for  pn  leaves 

ph(Vh) = 

If  the  above  equation  is  substituted  in  Eq.[2.116],  then 


Substituting  Eq.[2.118]  into  Eq.[2.115]  results  in 


(2.116) 


(2.117) 


(2.118) 


p(e,V)=p„(V)  1 


-TO 


r 

+  ve 


(2.119) 


The  reference  Hugoniot  state  including  the  shock  speed  has  4  unknowns,  i.  e.  the  speed  of  the  discontinuity 
U8,  vh,  &H)  and  pn •  There  are  only  three  Rankine-Hugoniot  equations,  thus  some  other  equation  must  be 
supplied  to  complete  the  system  of  unknowns  at  the  reference  Hugoniot2.  The  Experimental  Hugoniot  will 
be  used  for  this  purpose.  Combining  Eq.[2.37]  and  either  of  Eqs.[2.39]  or  Eq.[2.40]  results  in 

2 If  the  Mie-Grueneisen  were  designed  such  that  the  reference  pressure  were  a  function  of  energy  as  well  as  specific  volume, 
then  Eq.[2.116]-Eq.[2.118]  would  be  sufficient  to  define  the  reference  Hugoniot  state. 
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(2.120) 


Ph(V)  = 


4(Vo  -  V) 
[s(Vb  -  V)  -  E0]2 


Substituting  this  expression  into  Eq.[2.119]  leaves 


p(e,V) 


cUVo  ~  V) 

[s(V0  -V)-  Vo}2 


(2.121) 


which  is  the  Mie-Grueneisen  Equation  of  State.  This  empirical  condition,  where  s  and  cq  are  material 
constants  ,  is  particularly  useful  both  because  of  its  simplicity  and  that  it  has  been  experimentally  verified 
to  represent  many  materials  quite  well. 


2.2  Numerical  Aspects 

This  section  presents  the  numerical  algorithms  used  to  simulate  the  impact  dynamics  of  uniaxial  strain.  Sub¬ 
section  2.2.1  presents  the  methods  used  to  integrate  the  mass  continuity,  momentum,  and  energy  equations 
cast  in  an  Eulerian  reference  frame  using  a  strictly  conservative  integration  method.  Section  2.2.2  describes 
upwind  approaches  for  evaluating  the  flux  used  in  the  integration  of  the  conservative  equations  in  a  manner 
which  is  essentially  non-oscillatory  without  the  need  for  explicit  dissipative  dampening.  Section  2.2.3  derives 
the  eigenvalues  required  to  construct  the  characteristic  velocity  used  in  the  integration  scheme.  The  lasts 
section  2.2.5  describes  the  technique  developed  to  avoid  applying  the  traditional  “Turning  Method”  when 
the  numerical  algorithm  transitions  from  the  elastic  to  the  plastic  portion  of  the  surface. 

2.2.1  Integration  Scheme  for  the  Conservation  Equations 

The  equations  of  motion  as  presented  from  an  Eulerian  perspective,  do  not  implicitly  solve  for  the  deformation 
of  the  continuum,  but  rather  solve  for  the  unknown  velocity  field.  In  contrast,  Lagrangian  techniques  assume 
a  material  perspective  and  therefore  track  the  material  motion  as  the  deformation  occurs.  The  disadvantage 
of  the  latter  is  that  the  discretization  of  the  deformed  media  can  become  highly  skewed  resulting  in  reduced 
accuracy.  On  the  other  hand,  the  Eulerian  technique  allows  mass  to  pass  through  the  internal  cell  structures 
providing  the  opportunity  for  the  numerical  technique  to  adjust  the  nodal  placement  within  the  continuum. 
The  current  philosophy  is  to  require  that  nodes  which  lie  on  material  boundaries  remain  on  the  boundaries 
for  all  time,  whereas  the  interior  nodes  may  deform  in  any  manner  to  maintain  the  accurate  solution  of  the 
problem.  The  requirement  that  boundary  nodes  remain  on  the  boundary  enters  the  solution  as  a  boundary 
condition  where 


=  g{vi)  for  X{  €  [boundary]  (2.122) 

Note  that  g(v{)  for  the  multidimensional  case  does  not  require  the  boundary  nodes  to  move  at  local  particle 
velocities  (though  this  would  be  a  solution),  but  simply  requires  that  the  boundary  nodes  lie  on  the  surface 
of  the  deformed  material.  This  condition  is  also  used  to  specify  collisions  with  other  surfaces,  or  a  free 
boundary  which  evolves  with  the  solution.  Grid  motion  on  the  interior  given  by 
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dx  * 

— —  =  h(vi)  for  X{  €  [interior] 


(2.123) 


can  support  any  function  h(vi)  which  is  compatible  with  the  boundary  deformation  to  achieve  a  desirable 
nodal  distribution.  However  the  grid  advection  speed  contributes  to  the  eigenvalues  of  the  system  and  since 
the  maximum  time  step  is  inversely  proportional  to  the  maximum  eigenvalue  of  the  system,  functions  which 
move  the  interior  nodes  against  the  particle  motion  can  significantly  reduce  the  time  step  taken  by  the 
integration  algorithm.  For  this  reason,  care  should  be  taken  when  choosing  h(vi)  to  minimize  this  influence. 

While  the  Eulerian  form  of  the  equations  of  motion  provides  additional  flexibility  in  the  discretation  of 
the  domain,  this  flexibility  adds  complexity  to  the  integration  of  the  equations  of  motion.  To  illustrate  this, 
consider  the  divergence  form  of  the  one-dimensional  equations  of  motion  given  by 


dU  df 

dt  dx 


=  0 


The  equation  is  integrated  over  an  arbitrary  control  volume  resulting  in 


/,W 


df 

dx 


dV  =  0 


Gauss’s  Theorem  is  used  to  represent  the  second  term  as  a  surface  integral  giving 


(2.124) 


(2.125) 


f  ^dV  +  £(f  ■  rj)dS  =  0  (2.126) 

Here,  the  surface  S  bounds  the  control  volume  V,  and  fj  is  a  unit  outward  normal  on  the  surface  S.  The 
volume  and  the  control  surface  are  assumed  to  be  constant  during  integration,  therefore  Eq.[2.126]  can  not 
be  used  to  integrate  the  equations  of  motion  if  the  grid  is  deforming.  This  concept  and  the  issues  surrounding 
a  moving  grid  for  Eulerian  equations  lead  to  the  definition  of  the  so-called  Geometric  Conservation  Law[7] 
or  GCL.  A  recent  investigation  conducted  by  [8]  determined  that  if  the  the  numerical  solution  could  exactly 
reproduce  a  constant  flow  field  subjected  to  grid  motion  ,  then  the  solution  was  at  least  first-order  accurate 
in  time.  However,  no  options  were  offered  for  schemes  requiring  a  higher  level  of  temporal  accuracy.  Is  it 
therefore  possible  to  have  a  method  which  is  at  least  second-order  accurate  in  time  with  a  moving  mesh? 
Fortunately,  the  answer  is  that  any  order  of  temporal  accuracy  is  possible  if  the  integration  of  the  equations 
of  motion  properly  account  for  the  deformation  of  the  grid. 

The  fundamental  problem  of  integrating  Eq.[2.126]  subjected  to  a  moving  mesh  comes  as  the  result  of 
the  conventional  wisdom  of  separating  the  time  derivative  from  the  spatial  derivatives  before  the  application 
of  Gauss’  Theorem.  As  early  as  1968,  [9]  and  later  [10]  recognized  that  the  divergence  form  of  the  equations 
of  motion  could  include  time  as  a  coordinate  direction,  thus  the  nabla  or  del  operator  V  could  be  written  in 
the  three-dimensional  form  for  cartesian  components  as 

v  =  (Jh  +  Tx*x  +  +  I~ziz)  (2-127) 
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For  the  one-dimensional  problem,  the  volume  represents  an  (t,x)  area.  For  a  two-dimensional  problem  the 
volume  is  a  (t,x,y)  time-space  slab,  and  a  three-dimensional  problem  produces  a  (t,x,y,z)  four-dimensional 
volume.  With  this  operator,  the  conservation  equations  can  now  be  represented  in  the  simple  divergence 
form  as 


V*F  =  £([/)  (2.128) 

where  F  =  (t/,  f).  The  volume  integral  can  now  be  applied  over  Eq.[2.128]  with  deforming  geometry  without 
any  special  consideration  as  this  new  definition  of  volume  V  is  constant  over  the  integral  giving 


<f  (V  *  F)dV  =  (f  S(U)dV  =  S0V  (2.129) 

Jv  Jv 

To  evaluate  the  right-hand  side  of  Eq.[2.128],  the  well  known  Mean- Value  Theorem  is  applied  where  So  is 
the  value  of  S(U)  evaluated  at  some  interior  point  in  the  volume.  In  practice,  it  might  be  better  to  use  some 
form  of  numerical  quadrature  to  produce  a  more  accurate  result.  The  simplification  for  the  current  research 
is  better  justified  in  that  S(U)  =  0  in  the  absence  of  body  forces  which  are  neglected  in  this  study.  For 
completeness  however,  the  term  is  left  in  the  derivation. 

Applying  Gauss’s  Theorem  to  the  left-hand  side,  Eq.[2.129]  is  rewritten  for  the  one-dimensional  equation 
set  as 


/  (F  •  i ))dS  =  (f  (Urjt  +  ffjx)dS  =  50V 

Js  Js 


(2.130) 


where  fj  =  (^,  fjx)  is  a  unit  outward  normal  to  the  control  volume.  Applying  the  discrete  analogy  of  Eq. [2. 130] 
to  each  edge  of  the  control  volume  leaves 


Ew, + fjvxj)  iSji  =  J2-Fi  =  s°v 

j= i  i= i 

Expanding  Eq.[2.131]  on  each  surface  and  solving  for  the  term  with  [/3  gives 


(2.131) 


(U3  Ax)n+1  =  (Ui  Ax)n  +  T%-F2x  +  S$V  =  (U1Ax)n  +  At/"3  (2.132) 


where 


A U?}3  =  F2  -  T%  +  S^V  (2.133) 

Eqs.[2.132,2.133]  represent  a  first-order  O(At)  accurate  Euler  integration.  To  extend  the  method  to  second- 
order  accuracy,  a  corrector  step  is  added  which  is  a  temporal  adaption  of  the  well  known  MacCormack 
explicit  algorithm  [11]  which  is  a  member  of  the  Lax-Wendroff  method.  It’s  advantage  over  the  classic  RK 
class  of  schemes  is  that  the  predictor  step  calculates  to  the  n- 1-1  time  level,  not  at  some  fractional  increment. 
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As  a  result,  boundary  conditions  applied  after  the  predictor  but  before  the  corrector  step  are  straight  forward 
to  implement.  The  two  step  algorithm  is  as  follows; 

(UAxy  =  (UiAx)n  +  AU?  (2.134) 

(UiAx)n+1  =  (UiAx)*  +  i(AE/;  -  Atf,n)  (2.135) 

The  integration  algorithm  given  above  is  very  efficient  requiring  only  two  levels  of  memory;  one  level 
is  used  to  store  variables  at  time  level  n,  and  the  second  is  used  to  store  AU.  The  efficiency  is  gained  by 
storing  both  levels  of  A U  in  the  same  array.  After  U*  is  calculated  during  the  predictor  step,  the  array 
that  holds  A Un  is  negated.  On  the  corrector  step,  the  new  values  of  AU  are  simply  summed  into  the  array 
holding  A Un.  The  corrector  step  only  uses  U*  and  At/’s  therefore  U*  calculated  during  the  predictor  step 

can  overwrite  Un .  This  efficient  memory  overlap  scheme  has  the  added  benefit  that  the  same  routine  that 

calculates  A Un  will  also  calculate  AU*. 

The  linear  stability  requirement  for  this  algorithm  is 

ftjcjmax  <  CpL  <  j  (2.136) 

Ax 

where  |c|max  is  the  maximum  eigenvalue  of  the  system,  and  CFL  is  the  well  known  Courant-Friedrichs  and 
Lewy  number. 

2.2.2  Flux  Evaluation 

Most  upwind  schemes  vary  only  in  the  method  used  to  compute  the  interface  flux  such  as  J~2  and  in 
Eq.[2.133].  The  schemes  are  designed  to  damp  oscillations  in  the  vicinity  of  discontinuities  while  maintaining 
a  sharp  profile.  For  the  current  application,  the  dampening  of  oscillations  is  particularly  critical  for  the 
velocity  profile  where  the  velocity  behind  the  shock  is  zero.  A  small  amount  of  negative  velocity  will  cause 
the  stress  to  unload  thereby  creating  large  perturbations  in  the  stress  profile. 

Though  much  work  has  been  done  on  the  development  of  upwind  schemes  in  the  last  thirty  years  (See  [12] 
[13]  for  an  introduction),  most  of  these  can  be  characterized  as  using  one  of  two  fundamental  approaches  to 
flux  construction;  flux  interpolation  [14]  [15]  [22]  [17],  and  variable  extrapolation  followed  by  flux  evaluation [18]. 
Some  of  the  algorithms  having  been  developed  for  applications  in  Computational  Fluid  Dynamics  depend 
on  properties  not  applicable  for  the  current  application.  Roe’s  scheme,  for  example,  depends  on  the  flux 
function  F(U)  being  homogeneous  of  order  one  in  I/,  namely,  AF  =  U  which  does  not  hold  for  the  current 
equation  of  state[2.121]. 

Several  of  the  schemes  referenced  in  the  previous  paragraph  were  applied  with  some  success  to  the  current 
application.  The  CUSP  scheme  of  [15]  provided  good  results  if  the  user  was  willing  to  tune  the  flux  limiter 
for  the  specific  test  case.  In  general,  most  of  the  schemes  were  either  too  dissipative  to  prevent  the  overshoot 
of  shocks,  or  did  not  provide  enough  dissipation  resulting  in  massive  oscillations.  The  scheme  providing  the 
best  overall  performance  was  the  WENO  scheme  [19]  which  adaptively  computes  the  best  stencil  to  maximize 
smoothness  of  the  solution. 

Using  the  notation  of  [20],  we  general  T<i  and  as  occurring  on  a  material  interface  defined  by 

Ti+i=?t+t;+1  (2.137) 
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where 


ft  =  WEN05(rt-2,rt-vrt,rt+vft+2) 

(2.138) 

fr+i  =  WEN05(Fr+3,jr-+2,Fr+l,jr- 

(2.139) 

(2.140) 

and 

m  =  Max{\T'{Ui-i)\,  \!F'(Ui)\i  |.F'({/i+i)|) 

(2.141) 

where  \T'{Ui)\  is  the  maximum  eigenvalue  (section  2.2.3).  The  WEN05(a,b,c, 

d ,  e)  operator  defines  three 

interpolated  values 

,  a  7b  11 c 

91  “  3_6'  +  _6“ 

(2.142) 

b  5c  d 
q2  ~  "6  +  '6+3 

(2.143) 

c  5d  e 

qS  —  3  +  ~q  g 

(2.144) 

and  three  smoothness  functions 

/Si  =  13(a  -  2b  +  c)2  +  3 (a  -  46  +  3c)2 

(2.145) 

IS2  =  13(6  —  2c  -f  d)2  +  3(d  —  b)2 

(2.146) 

IS3  =  13 (c  -  2d  +  e)2  +  3(3c  -4 d  +  e)2 

(2.147) 

The  WEN05(a,  6,  c,  d,  e)  operator  is  then  defined  as 

xxt  i~>  /  L  J  \  Uiqi  +  UJ2q2  +  U)3q3 

TT£iV05(a,  o,c,  a,  e)  = - 

(2.148) 

W\  +  Ct>2  T  ^3 

where  the  weights  ui  are  defined  as 

1 

Wt  "  (e  +  ISi)2 

(2.149) 

For  all  computations,  e  =  10  6  as  suggested  by  [20]. 

The  truncation  error  for  the  discretation  is  of  order  0( Ax5)  in  space  and  of  order  0(&t2)  in  time  from 
Eqs[2.134,  2.135]. 
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2.2.3  Eigenvalue  Calculation 

The  Weighted  ENO  numerical  procedure  outlined  in  section  2.2.2  has  the  advantage  that  only  a  single 
wavespeed  is  needed  (see  Eq.[2.141]),  unlike  other  upwind  schemes  which  require  the  similarity  transformation 
matrices  such  as  the  Rieman  scheme  and  others.  For  a  complex  equation  of  state,  the  calculations  can  be 
quite  tedious  particularly  in  multiple  dimensions. 

The  wave  speed  used  in  Eq.[2.141]  is  the  maximum  eigenvalue  of  the  derivative  of  the  flux  T{U)  Eq.[2.130]. 
The  derivative  of  the  flux  is  given  by 


9T{U) 

dU 


q  f 

IVt  +  gjjVx 


or  in  matrix  form  as 


where 


dT{U) 

dU 


Vt 

(i?  -  *•')')* 


T)x  0 

Vt  +  uVx(2  -  ir0p0)  ^VxT'oPo 

^(P— S.+pet-  STopo}  rk+VxUT ^ 


(2.150) 


(2.151) 


dp  _  Cgpg[r0(/j  -  po)  +  spp  -  (1  +  a)]  _  uiy0 
dp  ~  \p{s  - 1)  -  sp0]3  P 

For  completeness,  the  following  derivatives  were  substituted  out  of  Eq. [2. 151] 

dp  _  -tdy0 

d(pu)  ~  p 

dp  _  TqPq 
dipet )  _  P 

The  eigenvalues  for  Eq.[2.151]  are 


Ai  =  vt  +  uVxA{^  =  Ai  ±cj]x 


where  c  is  defined  as 


(2.152) 


(2.153) 


(2.154) 


(2.155) 


c 


\p  ~  Sx  +  pjet  ~  M2)]ro/Jo 
P2 


(2.156) 


Note  that  this  expression  is  not  equivalent  to  the  sound  speed  associated  with  the  propagation  of  acoustic 
perturbations  through  the  material.  To  calculate  the  true  sound  speed,  it  is  necessary  to  represent  the 
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stress  and  strain  equations  as  partial  differential  equations  and  include  their  contributions  in  the  eigenvalue 
computation. 

To  complete  the  definition  of  defined  in  Eq.[2.141],  we  define 


=  l^itl  +  ci\Vxi\  (2.157) 

Note  that  if  the  grid  deforms  at  the  particle  velocity  rate  =  u),  then  r)t  =  —dx  =  - uSt ,  and  r]x  =  St 
for  an  outward  normal.  For  this  case,  Ai  =  0  thus  the  grid  motion  was  able  to  eliminate  the  advective 
component  of  the  flux  leaving  only  the  pressure  terms.  This  reduces  the  maximum  eigenvalue  which  in 
turn  reduces  the  amount  of  dissipation  needed  for  the  upwind  scheme  (Eq.2.140),  as  well  as  increasing  the 
overall  computational  time  step.  The  computational  advantage  is  not  dramatic  for  cases  where  u  <  c,  but 
does  validate  the  intuitive  notion  that  the  grid  should  as  a  general  rule  move  at  particle  velocity  rates  when 
possible. 


2.2.4  Integration  of  the  Constitutive  Equations 

The  numerical  integration  technique  for  the  conservation  equations  provides  the  mechanism  to  account  for 
grid  deformation  both  in  time  and  space.  To  accommodate  grid  motion  for  the  constitutive  equations 
[2.19d]  and  [2.19e],  a  coordinate  transformation  will  be  considered  where  t  =  t(r),x  =  ar(r,  C)  and  inversely 
r  =  r(<),  C  =  Applying  the  chain  rule  and  expanding  the  transformed  equations  leaves 


dSx 

dr 

dl 

dr 


du 

dl 


where  the  wave  speed  Ct  is  given  by 


(2.158) 

(2.159) 


Ct  = 


u  —  xT 
xc 


(2.160) 


Note  that  the  substitution  of  ^  in  Eqs.[2.158]  and  [2.159]  reproduces  the  original  equations  [2.19d] 

and  [2.19e]  except  that  t  in  the  original  equation  is  replaced  with  r.  Therefore  from  the  perspective  of 
numerically  solving  the  total  differential  equations,  the  transformation  may  be  ignored. 

To  integrate  Eqs.[2.19d,2.19e],  assume  the  equations  can  be  written  in  the  general  form 


D<f> 

Dt 


f(4>,  t ) 


(2.161) 


where  the  partial  derivative  on  the  right-hand  side  has  been  replaced  by  some  numerical  approximation  and 
the  function  f(<f>,  t)  may  or  may  not  be  an  explicit  function  of  <j)  and  t.  A  second  order  accurate  two-step 
numerical  solution  to  Eq.[2.161]  is  given  by 
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(2.162) 

(2.163) 


<!>*  =  </>n  +  WV) 

<t>n+ 1  =  r  +  \st{f(r,t+8t)-f(r,t)} 

Note  that  the  integration  procedure  is  identical  to  that  used  for  the  conservation  equations.  This  is  important 
in  that  (i)  the  dependent  variables  will  be  at  the  same  point  in  time  after  the  predictor  step  thereby  simplifying 
boundary  conditions,  and  (ii)  that  the  schemes  share  a  common  linear  stability  analysis. 

2.2.5  A  Turning  Algorithm 

The  classical  numerical  procedure  for  calculating  stress  and  stain  rate  assumes  that  the  material  response 
during  loading  is  elastic  until  the  effective  stress  reaches  the  yield  surface.  However,  the  numerical  procedure 
rarely  produces  a  stress  exactly  on  the  yield  surface  and  instead  generates  a  forbidden  stress  above  the 
yield  surface.  The  so-called  “Return  Methods”  treat  this  problem  by  geometrically  projecting  this  forbidden 
stress  onto  the  yield  surface  using  techniques  such  as  the  Radial  or  Oblique  Return  Method[21].  While  these 
methods  have  been  used  with  great  practical  success,  their  derivation  does  not  strictly  satisfy  the  differential 
equations  for  stress  and  strain  rate. 

The  conventional  interpretation  of  the  constitutive  equations  which  govern  the  stress  and  strain  pro¬ 
duction  is  that  the  transition  from  elastic  to  plastic  deformation  occurs  at  the  yield  surface  with  a  slope 
discontinuity.  This  mathematical  interpretation  is  commensurate  with  the  concept  that  the  material  follows 
either  an  elastic  or  plastic  response.  While  this  behavior  is  approximately  correct,  a  more  realistic  interpre¬ 
tation  is  that  the  transition  from  elastic  to  plastic  deformation  behaves  more  as  a  continuous  function  with 
a  pronounced  change  of  slope  at  the  yield  surface.  Based  on  this  interpretation,  a  new  method  is  presented 
which  assumes  that  the  material  behavior  can  be  partially  plastic  at  a  point  before  the  stress  violates  the 
yield  surface.  By  selecting  the  appropriate  level  of  plasticity  for  the  next  iteration,  the  algorithm  will  com¬ 
pute  the  value  of  stress  that  falls  exactly  on  the  yield  surface.  Since  a  forbidden  stress  is  never  used  in  the 
calculation,  this  method  is  coined  a  “Turning  Method”  as  it  provides  an  incremental  transition  from  the 
elastic  to  the  plastic  region3. 

If  the  current  effective  stress  is  below  the  yield  curve  and  in  loading  (point  E  as  shown  in  Figure. 2. 5) 
,  then  an  elastic  response  is  assumed  and  a  first  order  estimate  of  the  stress  is  calculated.  If  the  resulting 
effective  stress  exceeds  the  yield  surface,  then  the  first  order  estimate  is  discarded.  From  this  point,  the 
material  response  is  assumed  to  be  in  a  mixed  elastic/plastic  mode.  To  update  to  the  new  value  of  stress 
and  strain,  the  “Turning  Method”  utilizes  the  existing  stress  and  strain  rate  equations  [2.19d]  and  [2.19e], 
except  that  the  parameter  k  defined  in  Eq.[2.10]  is  allowed  to  take  on  a  fractional  value.  Casting  Eq.[2.19d] 
and  Eq.[2.16]  in  an  incremental  form  and  substituting  Eq.[2.10],  we  have 


where 


1  /»Ac 7 

2  GuxAt 


(2.164) 


Act  —  H(ttrial)  O elastic  (2.165) 

3 An  early  temptation  was  to  refer  to  this  scheme  as  the  “Method  of  No-Return”.  Fortunately,  good  taste  and  a  more 
descriptive  name  won  out. 
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and 


/.  =  Sgn(Sx)  (2.166) 

where  the  value  k  is  assumed  constant  during  this  elastic/plastic  phase.  If  H'(e)  is  zero,  then  the  response  is 
purely  plastic  and  A <7  is  constant  resulting  in  a  non-iterative  scheme.  If  H'(e)  is  non-zero,  then  etriai  is  set 
to  zero  in  Eq.[2.165],  A  a  is  calculated,  followed  by  k  Eq.[2.164],  then  fi  Eq.[2.10].  The  process  is  repeated 
until  the  change  in  k  is  sufficiently  small  which  in  practice  requires  only  a  few  iterations. 

Another  way  to  achieve  the  same  result  is  to  consider  that  the  material  is  plastic  (namely,  A;  =  1),  but 
that  the  slope  of  the  material  H'  can  be  increased  beyond  what  is  physical  to  transition  to  the  yield  surface. 
This  equivalent  slope  can  be  calculated  using  /3  from  the  current  scheme  and  solving  Eq.[2.10]  for  H' .  This 
illustrates  that  the  scheme  simply  clips  the  edge  of  the  stress-strain  curve  in  order  to  avoid  violating  the 
yield  surface. 


2.3  Numerical  Results 

To  illustrate  the  accuracy  of  the  numerical  algorithm,  the  test  cases  presented  in  the  examples  will  be  solved 
numerically  and  compared  to  the  exact  solutions.  The  Figures  for  the  examples  are  shown  in  Chapter  6. 

In  order  to  compute  the  numerical  solution  for  the  hydrostatic  and  elastic  solution  using  the  experimental 
hugoniot,  it  was  necessary  to  find  a  function  which  expresses  pressure  as  a  function  of  known  properties, 
which  in  effect,  provides  the  role  of  an  equation  of  state.  To  find  a  suitable  functin,  consider  the  expression 
for  a  Raleigh  line  given  by 


%—^  =  -pl{Us-v2f  (2.167) 

which  connects  some  initial  point  1  on  a  p  —  V  Hugoniot  to  any  other  point  2  on  the  same  p  —  V  curve. 
Using  this  function,  a  plot  of  pressure  versus  particle  speed  can  be  produced  that  shows  all  possible  states 
after  passing  through  a  single  shock.  Substituting  the  equation  for  the  right-running  experimental  hugoniot 
and  the  jump  condition  for  continuity,  and  after  considerable  manipulation  results  in 


A  _  PlpTClipi  -p2) 
K/>l-p2)  +  />2]2 


(2.168) 


For  the  hydrostatic  case,  deviatoric  stress  will  be  0.  Assuming  the  reference  pressure  is  also  0  and  that  the 
reference  density  is  p0y  results  in 


p_  PopC$(p-Po) 
[«(a>  -  p)  +  p]2 


or  expressed  in  terms  of  specific  density  as 


CjjVo  -  V) 

[V0(s  -  1)  -  sV}2 


(2.169) 


(2.170) 
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For  the  elastic  case  using  the  experimental  hugoniot,  the  effects  of  the  deviatoric  stress  must  be  included. 
Assuming  the  reference  deviatoric  stress  is  0  leaves 


P  _  a  PoPCqCp-Po) 

r-^  +  [s(p0-p)+p}l 


or  expressed  in  terms  of  specific  density  is 


(2.171) 


P  =  SX  + 


cgw,  -  v) 

[Vb(s  -  1)  -  sV]2 


(2.172) 


As  described  in  the  examples,  the  configuration  of  the  test  cases  are  the  impact  of  two  metal  bars  with 
properties  given  in  Table  2.1.  Their  initial  pressure,  energy,  and  deviatoric  stress  are  zero,  i.  e.  pi  =  e i  = 
SXl  =  0.0  and  p4  =  e 4  =  SXA  =  0.0  in  both  bars,  and  their  initial  densities  are  8930.0  Kg/m3.  The  velocities 
of  the  bars  are  given  in  the  context  of  each  test  case. 


2.3.1  Hydrostatic  Test  Case 

See  the  example  in  Section  2.1.2  for  details  of  this  test  case.  The  results  of  velocity  vs.  pressure,  density,  and 
energy  are  shown  in  Figures  [6.1],  [6.2],  and  [6.3],  respectively.  Table  2.5  shows  a  comparison  of  the  results. 
Only  state  2  is  shown  since  state  3  is  identical  to  state  2. 

The  results  of  this  experiment  indicate  that  the  numerical  algorithm  is  doing  a  very  good  job  of  capturing 
the  fundamental  physics  of  the  problem. 


Table  2.5:  Exact  vs.  Computed  Results  for  Hydrostatic  case(State  2). 


Parameter  (Units) 

Computed 

Exact 

|  Error  | 

Velocity(  m/s) 

0.0 

0.000 

0.0% 

Pressure(  GPa) 

3.652 

3.653 

0.0% 

Density (  Kg/m3) 

9154.0 

9153.866 

0.0% 

Energy  (  J/Kg) 

4991.0 

5000.000 

0.2% 

Shock  Vel.(  m/s) 

-3994.0 

-3989.000 

0.1% 

2.3.2  Elastic  Solution  using  the  Experimental  Hugoniot. 

See  the  example  in  Section  2.1.2  for  details  of  this  test  case.  The  results  of  velocity  vs.  pressure,  density, 
energy,  and  deviatoric  stress  are  shown  in  Figures  [6.4],  [6.5],  [6.6],  and  [6.7],  respectively.  The  comparison 
between  the  exact  solutions  and  the  numerical  solution  are  shown  in  Table  2.6. 

This  test  case  illustrates  some  oscillation  particularly  with  the  energy  and  the  deviatoric  stress  near  the 
impact  point.  However,  the  scheme  is  able  to  capture  the  correct  values  away  from  the  oscillations. 


2.3.3  Elastic  Solution  using  the  Mie-Grueneisen  Equation  of  State 

See  the  example  in  Section  2.1.2  for  details  of  this  test  case.  The  results  of  velocity  vs.  pressure  ,  density, 
energy,  and  deviatoric  stress  are  shown  in  Figures  [6.8],  [6.9],  [6.10],  and  [6.11],  respectively.  The  comparison 
between  the  exact  solutions  and  the  numerical  solution  are  shown  in  Table  2.7. 
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Table  2.6:  Exact  vs.  Computed  Results  for  Elastic  case(State  2)  using  the  Experimental  Hugoniot. 


Parameter(Units) 

Computed 

Exact 

|  Error  | 

Velocity  (  m/s) 

0.0 

0.0 

0.0% 

Pressure(  MPa) 

200.09 

201.319 

0.6% 

Density(  Kg/m3) 

8952.65 

8952.637 

0.0% 

Energy  (  J/Kg) 

50.0 

49.9 

0.1% 

Deviatoric  Stress  (  MPa) 

151.5837 

-152.123 

0.3% 

Normal  Stress(  MPa) 

-353.1726 

-353.1726 

0.0% 

Shock  Vel.(  m/s) 

-3994.0 

-3967.8 

0.5% 

This  test  case  is  similar  to  the  previous  case  in  that  both  exhibit  oscillations  in  deviatoric  stress  and 
energy  near  the  impact  point.  However,  like  before  the  computed  values  are  very  close  to  the  exact  values 
away  from  the  oscillations. 


Table  2.7:  Exact  vs.  Computed  Results  for  Elastic  case(State  2)  using  the  Mie-Grueneisen  Equation  of 
State. _ 


Parameter(Units) 

Computed 

Exact 

|  Error) 

Velocity  (m/s) 

0.0 

0.0000 

0.0% 

Pressure  (MPa) 

295.4 

295.2880 

0.0% 

Density  (Kg/m3) 

8948.9 

8948.9250 

0.0% 

Energy  (J/Kg) 

49.9 

50.0000 

0.1% 

Deviatoric  Stress  (MPa) 

-127.0 

-126.9700 

0.0% 

Shock  Vel.  (m/s) 

-4741.4 

-4718.6411 

0.4% 

2.3.4  Two- Wave  Solution  Using  the  Mie-Grueneisen  Equation  of  State. 

See  the  example  given  in  Section  2.1.3  for  details  of  this  test  case.  The  results  of  velocity  vs.  pressure, 
density,  energy,  and  deviatoric  stress  are  shown  in  Figures  [6.12],  [6.13],  [6.14],  and  [6.15],  respectively.  The 
comparison  between  the  exact  solution  and  the  numerical  solution  is  shown  in  Table  2.8. 

This  test  case  has  a  much  higher  impact  velocity  than  the  previous  case  much  as  indicated  by  the  higher 
error,  although  the  oscillations  are  much  better  damped  in  proportion  to  the  previous  case. 


Table  2.8:  Exact  vs.  Computed  Results  for  the  Two- wave  Shock  using  the  Mie-Grueneison  Equation  of 
State. _ 


State  1' 

State  2 

Parameter  ( U  nits) 

Computed 

Exact 

Error 

Computed 

Exact 

|  Error  | 

Velocity  (m/s) 

84.502 

84.234 

0.3% 

0.0 

0.000 

0.0% 

Pressure  (MPa) 

457.360 

466.747 

2.0% 

3559.5 

3558.695 

0.0% 

Density  (Kg/m3) 

8959.270 

8959.826 

0.0% 

9147.8 

9147.919 

0.0% 

Energy  (J/Kg) 

120.369 

124.272 

3.1% 

5199.0 

5219.884 

1.9% 

Elastic  Shock  (m/s) 
Plastic  Shock  (m/s) 

-4678.3 

-4635.9538 

0.9% 

-3998.5 

-4012.5352 

0.4% 

44 


2.3.5  Single- Wave  Plastic  Solution  using  the  Mie-Grueneisen  Equation  of  State. 

See  the  example  given  in  Section  2.1.4  for  details  of  this  test  case.  The  results  of  velocity  vs.  pressure, 
density,  energy,  and  deviatoric  stress  as  shown  in  Figures  [6.16],  [6.17],  [6.18],  and  [6.19],  respectively.  The 
comparison  between  the  exact  solution  and  the  numerical  solution  is  shown  in  Table  2.9. 

This  solution  exhibits  some  oscillation  as  the  result  of  a  very  small  undershoot  of  the  velocity  behind 
the  shock.  The  undershoot  causes  the  bar  to  unload  elastically,  then  an  oscillation  sets  up  when  the  stress 
attempts  to  recover  to  the  yield  surface.  One  possible  way  to  correct  this  behavior  is  to  use  a  higher  order 
temporal  algorithm  for  the  integration  scheme. 


Table  2.9:  Exact  vs.  Computed  Results  for  the  Plastic  Single- Wave  response  using  the  Mie-Grueneisen 
Equation  of  State. _ 


Parameter(Units) 

Computed 

Exact  | 

Error  | 

Velocity  (m/s) 

0.0 

0.000 

0.0% 

Pressure  ( GPa ) 

25.8 

25.800 

0.0% 

Density  (Kg/m3) 

10190.0 

10189.950 

0.0% 

Energy  (J/Kg) 

180073.3 

180458.200 

0.2% 

Deviatoric  Stress  (MPa) 

-200.0 

-200.000 

0.0% 

Shock  Vel.  (m/s) 

-4249.0 

-4135.954 

2.0% 
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Elastic  Response(k=0) 


Figure  2.5:  Figure  illustrating  the  “Turning  Method1 


Chapter  3 

Two-Dimensional  Considerations 


This  chapter  will  present  the  equations  of  motion  and  the  development  of  numerical  algorithms  to  solve  the 
equations  of  motion  for  hydrodynamic  simulations.  Two  examples  will  be  shown  demonstrating  the  use  of 
the  algorithms. 


3.1  Conservation  Equations 

In  the  absence  of  heat  conduction  and  any  sources  the  conservation  equations  for  two-dimensional  plane 
stress  is  written  in  the  vector  form  as 


8U  OF  8G 
dt  +  dx  +  dy 


(3.1) 


or  in  the  divergence  form  of 


V  •  F  =  0  (3.2) 

where 


F  =  (U,F,G) 


(3.3) 


and 
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pu 

pv 

pu 
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pu2  +  p-  Sx 
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PUV  -  Sxy 
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pUV  -  SX  y 

PV2  +P~  Sy 

per 

_  puer  +  u(p  -  Sx)  -  vSxy 

_  puer  +  v(p  -  Sy)  -  vSxy 

(3.4) 
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3.2  Constitutive  Equations 

The  constitutive  equations  written  to  satisfy  material  frame  indifference  for  plane  strain  is 

q  ( du  o/'~*  (  2  du  1  dv  3  fj  5x  i  \ 

b*y  (dt  ~  m  )  +  ZU  [  3  Qy  ~  2  PgSVxyl 

where 


I? 

£>£ 


e 


4>. 


xy 


-S  ^  +  s  —  +  s 

~bxdx+bydy+bxy 


(3.6) 


The  conservation  equations  (Section  3.1)  solved  with  the  constitutive  equations  and  closed  with  an  equation 
of  state  (Section  2.1.5)  form  the  complete  equations  of  motion  for  hydrodynamics. 


3.3  Integration  of  the  Conservation  Equations 

The  method  of  Integration  follows  closely  with  that  developed  in  Section  2.2  for  the  case  of  uniaxial  strain. 
Assuming  two  space  dimensions  Eq.[2.130]  is  written  as 


4  Gf}y)dS  —  0 


Expanding  the  terms  in  the  expression  above  we  have 


(3.7) 


^(U,7))  =  F-f}  = 


pUn 

puUn  4  (p  ~~  Sx)fjx  -  SXyfjy 
pvUn  ~  SXyfjx  +  (P  “  Sy)ffy 
peTUn  4  [u{p  —  Sx)  ~~  vSXy]f}x  4  [v{p  ~  sy)  ~  ^^xylVy 


(3.8) 


where  Un  =  {fjt  4  ufjx  4  vfjy). 

The  continuum  is  discretized  using  an  unstructured  grid  system  where  the  field  parameters  are  stored  at 
the  vertices  of  the  grid.  A  control  volume  (c.v.)  is  constructed  about  each  vertex  using  the  dual  medium 
grid  as  shown  in  Figure  [3.1]. 

Summing  Eq.[3.8]  on  all  faces  and  setting  the  total  flux  of  the  time-space  slab  to  zero  requires  that 


na 

Un+ Mn+1  -  UnAn  4  =  0 

i- 1 


(3.9) 
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Figure  3.1:  The  control  volume  using  a  dual  medium  grid. 


where  ns  is  the  number  of  lateral  sides  on  the  control  volume.  For  each  edge  connecting  the  center  node  “o”, 
and  each  neighboring  vertex,  two  flux  evaluations  are  required  corresponding  to  the  two  segments  connecting 
the  center  of  the  adjacent  cells  “c”  to  the  center  of  the  edge  “m”  Although  other  control  volumes  can  be 
constructed  such  as  connecting  the  center  of  the  cells  directly,  the  dual  medium  has  been  shown  to  provide 
the  best  results,  particularly  when  the  cell  structure  is  skewed  [22]. 

A  second  order  scheme  in  time  and  space  is  formed  by  a  two-step  explicit  integration  algorithm.  The 
predictor  step  is  a  simple  Euler  integration  derived  by  solving  Eq.[3.9]  for  the  term  at  the  “n  +  1”  level  or 


U*A*  =  UnAn  +  MJn 


(3.10) 


where 


A£/"  =  -XVf(U,^| 

i=l 


(3.11) 


The  corrector  phase  is  calculated  as 


un+lAn+l  _  vnAn  +  I[A jj*  _  Ajyn] 


(3.12) 


where 
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(3.13) 


ait  =  -]T;f;(u,7))|7?| 

i=l 

The  lateral  flux  component  T{U,f})  is  needed  at  the  middle  of  the  edge  segment  “m”,  yet  the  field 
parameters  are  stored  at  the  vertices.  There  are  many  methods  which  can  be  used  to  estimate  the  flux  at 
the  center  of  the  edge,  but  one  of  the  most  popular  is  the  class  of  upwind  methods.  Reference  [22]  suggests 
a  method  based  on  an  approximate  Rieman  method  where 

HVMo  =  E\\V\ (V(Uo,>7o)  +  - u0))  (3.14) 

The  absolute  value  sign  around  the  derivative  indicates  some  form  of  eigenvalue  decomposition  and  nor¬ 
malization.  This  procedure  requires  not  only  the  eigenvalue  calculation,  but  also  similarity  transformation 
matrices.  However,  the  complexity  of  this  operation  is  not  warranted.  The  sound  speed  is  much  larger  than 
the  local  particle  velocity  and  is  therefore  dominated  by  the  u  +  c  eigenvalue.  As  a  result,  the  absolute 
value  operation  in  Eq.[3.14]  is  replaced  with  a  scalar  velocity  a  which  satisfies  a  positivity  condition  on  the 
eigenvalues  where 

*>\\\i\max  (3-15) 

where 

A i  :  eigenvalues  of  ^  (3.16) 

m 

The  eigenvalue  calculation  is  given  in  Section  [2.2.3].  Substituting  a  in  Eq.[3.14]  leaves 

•F(U,7))|0  =  fj  \\n\  (V(U0,7)o)  +  -  a(U  i  -  U0))  (3.17) 

3.4  Integration  of  the  Constitutive  Equations 

The  fundamental  concept  of  integrating  the  constitutive  equations  [3.5]  is  to  replace  the  partial  derivatives 
with  finite  volume  expressions,  then  applying  some  standard  technique  for  integrating  the  resulting  first- 
order  Ordinary  Differential  Equation.  The  integration  technique  for  two-dimensions  is  the  same  as  shown  in 
Section  [3.2]  and  will  not  be  repeated  here. 

3.5  Results  for  Two-Dimensional  Plane  Strain 

The  concepts  explored  in  the  case  of  uniaxial  strain  were  extended  into  two-dimensions  primarily  to  under¬ 
stand  if  the  numerical  algorithms  would  extend  into  multi-dimensions  as  expected.  Two  simple  test  cases 
are  presented  that  illustrate  reasonable  qualitative  plastic  behavior,  though  the  results  were  not  validated. 
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The  first  test  case  was  a  finite  length  bar  impacting  a  non-deformable  surface  at  1000  m/s  (See  Figure 
[6.20]).  The  boundary  condition  on  the  left  wall  was  zero  horizontal  deflection,  thereby  simulating  a  cen¬ 
terline.  The  boundary  condition  imposed  on  the  bottom  surface  was  no  vertical  deflection.  The  right  and 
top  boundaries  were  free  to  deform  as  calculated  and  imposed  a  zero  stress  boundary.  Ideally,  the  boundary 
condition  should  be  zero  surface  traction,  i.  e.,  njcTji  =  0  but  the  details  of  the  implementation  were  never 
explored.  In  addition,  the  nodes  on  the  interior  of  the  simulation  were  moved  during  each  time  step  such 
that  LaPlace’s  solution  was  approximately  satisfied  at  each  time  step,  or 


V2z  +  V2y  =  0  (3.18) 

This  adaptivity  provides  a  smooth  variation  of  the  interior  grid  nodes  to  large  deformations  at  the  boundary. 
Similar  results  were  obtained  if  the  interior  nodes  were  allowed  to  follow  the  particle  velocity  rate. 

The  second  test  case  is  identical  to  the  first  test  case  except  a  notch  was  included  on  the  output  boundary. 
The  initial  configuration  is  shown  in  Figure  [6.21],  and  the  deformed  geometry  is  shown  magnified  in  Figure 
[6.22]. 

Neither  of  these  two  test  cases  were  evaluated  in  any  detail  other  than  to  observe  that  the  results 
indicated  basic  plastic  response  (i.e.  mushrooming  of  the  base)  with  an  elastic  rebound  which  was  observed. 
The  boundary  conditions  were  not  accurate  and  therefore  would  never  compare  favorably  with  data. 
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Chapter  4 

Conclusion 


The  overall  conclusion  is  an  Eulerian  form  of  the  equations  of  motion  for  hydrodynamics  can  be  used  to 
numerically  simulate  a  wide  wide  of  impact  velocities  with  great  accuracy  with  little  oscillation  in  the 
solution.  The  best  observed  upwind  technique  to  solve  these  equations  is  the  ENO  technique  which  produces 
the  sharpest  shocks  with  the  least  amount  of  smearing.  However,  very  high  orders  of  accuracy  in  both  time 
and  space  are  necessary  to  achieve  accurate  solutions  without  oscillations,  although  correct  jump  conditions 
can  be  achieved  in  the  presence  of  oscillations  with  lower  order  solutions.  Small  errors  in  velocity,  particularly 
shock  undershoots,  can  cause  the  deviatoric  stress  to  unload  generating  large  errors  and  oscillations  in  the 
stress  field. 

Exact  solutions  to  impact  problems  using  multiple  materials  can  be  produced  using  the  Mie-Grueneisen 
Equation  of  State,  though  not  all  the  solutions  are  explicit  and  must  be  solved  using  numerical  techniques 
to  handle  algebraic  non-linearity. 

It  is  not  desirable  to  solve  the  conservation  equations  and  the  constitutive  equations  as  a  coupled  set 
of  partial  differential  equations.  It  was  illustrated  that  there  were  no  discovered  divergence  forms  that 
produced  eigenvectors  which  remain  bounded,  and  therefore  retained  its  hyperbolic  character.  Secondly, 
no  forms  could  be  found  which  produced  any  different  eigenvalues  beyond  that  of  the  conservative  system, 
therefore  no  gain  would  be  attributed  to  the  coupling.  Perhaps  the  strongest  argument  is  that  some  of 
the  exact  solutions  demonstrate  that  the  solution  to  the  conservative  equations  were  dependent  only  on  the 
normal  stress,  and  not  on  the  individual  stress  components1.  Therefore,  it  would  not  make  physical  sense  to 
couple  the  conservative  and  constitutive  equations  when  the  physics  are  not  coupled. 


1  As  an  example,  the  solution  for  the  hydrostatic  and  the  elastic  solution  produce  the  same  normal  stress  independent  of  the 
choice  of  stress  model 
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Chapter  5 

Recommendations 


The  research  effort  in  examining  the  potential  of  using  the  Eulerian  form  of  the  equations  of  motion  coupled 
with  an  upwind  solution  methodology  is  a  success.  The  result  of  the  numerics  illustrates  an  exceptional 
capability  to  capture  the  correct  jump  conditions  across  shocks  without  the  use  of  artificial  dissipation 
to  maintain  stability.  Although  shocks  are  reasonably  sharp  using  the  WENO  integration  algorithm,  the 
numerical  solution  for  very  high  impact  velocities  exhibits  a  tendancy  for  the  deviatoric  stress  to  unload  in 
the  presence  of  small  velocity  undershoots  downstream  of  the  shock.  This  defect  in  stress  in  turn  corrupts 
the  calculation  of  plastic  strain.  This  effect  may  be  minimized  by  using  a  scheme  with  fourth-order  or  higher 
accuracy  in  time  (the  current  scheme  is  fifth  order  accurate  in  space,  and  second  order  in  time). 

Another  strategy  which  needs  to  be  examined  is  the  level  of  improvement  in  the  solution  accuracy  that 
could  be  realized  by  applying  upwind  techniques  to  the  equations  of  motion  cast  in  the  Lagrange  form. 
Simple  numerical  experiments  using  central  difference  approximations  indicate  the  ability  to  capture  shocks 
sharply,  although  artificial  diffusion  has  to  be  applied  to  maintain  stability,  which  taints  the  accuracy  of  the 
solution.  Upwind  techniques  should  further  refine  the  shock  without  the  explicitly  added  artificial  viscosity 
resulting  in  higher  accuracy. 

Much  work  is  yet  to  be  done  on  both  the  numerical  model  and  the  validation  of  the  two-dimensional 
computer  code.  The  current  implementation  needs  to  be  expanded  to  include  the  axisymmetric  form  for 
validation  against  the  Taylor  impact  experimental  data. 
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Chapter  6 

Additional  Figures 
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Figure  6.1:  Pressure  for  hydrostatic  test  case. 
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Figure  6.2:  Density  for  hydrostatic  test  case. 
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Figure  6.3:  Energy  for  hydrostatic  test  case. 
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Figure  6.4:  Pressure  for  elastic  test  case  using  the  Experimental  Hugoniot. 
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Figure  6.6:  Energy  for  elastic  test  case  using  the  Experimental  Hugoniot. 
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Figure  6.7:  Deviatoric  Stress  for  elastic  test  case  using  the  Experimental  Hugoniot. 
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Figure  6.10:  Energy  for  elastic  test  case  using  the  Mie-Grueneisen  equation  of  state. 
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Figure  6.11:  Deviatoric  stress  for  elastic  test  case  using  the  Mie-Grueneisen  equation  of  state. 
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Figure  6.14:  Energy  for  two- wave  test  case  using  the  Mie-Grueneisen  equation  of  state. 
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Elastic/Plastic  Impact 


Figure  6.19:  Deviatoric  stress  for  single- wave  test  case  using  the  Mie-Grueneisen  equation  of  state. 
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Copper  Bar 
Impact  Vel.  1000  m/s 


Figure  6.20:  Pressure  for  impact  of  Copper  Bar  at  1000  m/s 
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Figure  6.21:  Initial  Configuration  of  notched  bar. 
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Figure  6.22:  Deformed  Configuration  of  notched  bar  with  Pressure. 
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